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ABSTRACT 

We present detailed clustering measurements from the 2dF QSO Redshift Survey 
(2QZ) in the redshift range 0.8 < z < 2.1. Using a flux limited sample of ~ 14,000 
objects with effective redshift Zcs = 1.47, we estimate the quasar projected correla- 
tion function for separations 1 < r//i~^Mpc < 20. We find that the 2-point correlation 
function in real space is well approximated by a power law with slope 7 = 1.5 ± 0.2 
and comoving correlation length tq = 4.8^15 ^~^ Mpc. Splitting the sample into three 
subsets based on redshift, we find evidence for an increase of the clustering amplitude 
with lookback time. For a fixed 7, evolution of rg is detected at the 3.6(7 confidence 
level. The ratio between the quasar correlation function and the mass autocorrelation 
function (derived adopting the concordance cosmological model) is found to be scale 
independent. For a linear mass-clustering amplitude CTg = 0.8, the "bias parameter" 
decreases from b ~ 3.9 at Zoff = 1-89 to 6 ~ 1.8 at Zcs — 1.06. From the observed 
abundance and clustering, we infer how quasars populate dark-matter haloes of dif- 
ferent masses. We find that 2QZ quasars sit in haloes with M > 10^'^Mq and that 
the characteristic mass of their host haloes is of the order of 10^^ Mq. The observed 
clustering is consistent with assuming that the locally observed correlation between 
black- hole mass and host-galaxy circular velocity is still valid at z > 1. From the 
fraction of haloes which contain active quasars, we infer that the characteristic quasar 
lifetime is iq ~ a few x 10^ yr at z ~ 1 and approaches 10® yr at higher redshifts. 

Key words: galaxies: active - galaxies: clustering - quasars: general - cosmology: 
theory - large-scale structure - cosmology: observations 



1 INTRODUCTION 

Recent dynamical studies have provided strong evidence for 
the existence of supermassive black holes in the centre of 
most nearby galaxies (for a review, see e.g. Richstone et al. 
1998). The mass of the central black hole seems to corre- 
late with the luminosity and the velocity dispersion of the 
spheroidal stellar component (e.g. Magorrian et al. 1998; 
Gebhardt et al. 2000; Ferrarese & Merritt 2000; Tremaine 
et al. 2002). The surprising tightness of the latter relation 
suggests the existence of a strong connection between the 
formation of supermassive black holes and the assembly of 
galactic spheroids (Silk & Rees 1998; Haehnelt & Kauff- 
mann 2000; Kauffmann & Haehnelt 2000; Monaco, Salucci 
& Danese 2000; Granato et al. 2004). 

The mounting evidence for the presence of supermassive 
black holes in nearby galaxies supports the theoretical belief 
that quasars are powered by black-hole accretion (Salpeter 
1964; Zel'dovich & Novikov 1964; Lynden-Bell 1969). For 
instance, the locally estimated mass density in black holes 



and the observed evolution of the quasar luminosity func- 
tion seem to be consistent with this hypothesis (Haehnelt, 
Natarajan, Rees 1998; Fabian & Iwasawa 1999; Salucci et al. 
1999; Yu & Tremaine 2002; Wyithe & Loeb 2003; Marconi 
et al. 2004). However, a detailed understanding of the phys- 
ical processes leading to quasar activity (and their connec- 
tion with galaxy formation) is still lacking. For this reason, 
even simple phenomenological models which are able to re- 
produce the observational results by selecting which cosmic 
structures could harbour quasars are of paramount impor- 
tance. 

In the currently favoured cosmological model, galaxies 
are expected to form within extended dark-matter haloes. 
At every epoch, the number density and clustering proper- 
ties of the haloes can be readily (and reliably) computed 
as a function of their mass. It is therefore of great interest 
to try to establish a connection between these haloes and 
different classes of cosmic objects. Even though the distri- 
bution of light sources within haloes is determined by com- 
plex physics, some of its properties can be computed with 
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a purely statistical approach. For instance, one can use the 
mean density and the clustering amplitude of a population 
of cosmic objects to determine the lowest-order moments of 
the "halo occupation distribution function", Pjv(M), which 
gives the probability of finding a given number of sources in 
a halo of mass M (e.g. Scoccimarro et al. 2001 and references 
therein). A number of "halo models" have been presented in 
the literature. These have been successfully used to describe 
the adundance and clustering properties of galaxies at both 
low (Peacock & Smith 2000; Seljak 2000; Scoccimarro et al. 
2001; Marinoni & Hudson 2002; Berlind & Weinberg 2001; 
Yang, Mo & van den Bosch 2003; van den Bosch, Yang & 
Mo 2003; Zehavi et al. 2003; Magliocchetti & Porciani 2003) 
and high (Bullock, Wechsler & Somerville 2002; Moustakas 
& Somerville 2002; Hamana et al. 2004; Zheng 2004) red- 
shift. One of the main goals of this paper is to use a similar 
approach to investigate how quasars populate dark matter 
haloes. As previosly discussed, this requires an accurate de- 
termination of the clustering properties of bright QSOs. 

Since the first detection of quasar clustering, (Shaver 
1984; Shanks et al. 1987), a number of surveys (continuously 
improved in terms of homogeneity, completeness and size) 
have been used to measure the 2-point correlation function 
of bright QSOs (lovino & Shaver 1988; Andreani & Cristiani 
1992; Mo & Fang 1993; Shanks & Boyle 1994; Andreani et 
al. 1994; Groom & Shanks 1996; La Franca, Andreani & 
Cristiani 1998; Grazian et al. 2004). The emerging picture 
is that quasars at z ~ 1.5 have a correlation length ro ~ 5 — 6 
h~^ Mpc, similar to that of present-day galaxies. It still is a 
matter of debate, however, whether ro significantly evolves 
with redshift (lovino & Shaver 1988; Groom & Shanks 1996; 
La Franca, Andreani & Cristiani 1998). This uncertainty is 
due to the joint effects of cosmic variance and small-number 
statistics: given the sparseness of the quasar distribution, a 
typical sample includes from a few-hundred to a thousand 
objects. In consequence, clustering is generally detected at 
a relatively low significance level (3 — 4 a). 

The development of efficient multi-object spectrographs 
has recently made possible a new generation of wide-area 
redshift surveys. Both the completed Two-degree Field 
(2dF) QSO Redshift Survey (Groom et al. 2004) and the 
on-going Sloan Digital Sky Survey (SDSS) Quasar Survey 
(Schneider et al. 2003) list redshifts for tens of thousands 
of optically selected quasars. A preliminary data release of 
the 2QZ has been used to estimate the evolution and the lu- 
minosity dependence of the quasar 2-point correlation func- 
tion in redshift-space (Groom et al. 2001; 2002). The final 
catalogue has been employed to measure the quasar power 
spectrum out to scales of 500 h~^ Mpc (Outram et al. 2003; 
see also Hoyle et al. 2002) and to constrain the cosmologi- 
cal constant from redshift-space distortions (Outram et al. 
2004). 

In this paper, we study the clustering properties of 
~ 14, 000 quasars extracted from the complete 2dF QSO 
Redshift Survey. In particular, we compute the evolution of 
their projected 2-point correlation function. This quantity 
is not affected by the distortion of the clustering pattern 
induced by peculiar motions as it measures the clustering 
strength as a function of quasar separation in the perpen- 
dicular direction to the line of sight. Using the largest quasar 
sample presently available, we are able to accurately mea- 
sure the real-space clustering amplitude in three redshift 



bins. Our results reveal a statistically significant evolution 
of the clustering length with redshift. 

We then combine our clustering study with the halo 
model to infer the mean number of optically selected quasars 
which are harboured by a virialized halo of given mass (the 
halo occupation number) and the characteristic quasar life- 
time. Our results can be directly used to test physical models 
for black-hole accretion. 

The layout of the paper is as follows. In Section |21 we 
present our quasar samples, we measure their projected cor- 
relation function and estimate the corresponding bias pa- 
rameters. In Section El we introduce the halo model and 
we discuss how the the halo occupation number, N{AI), 
is constrained by the observed abundance and clustering 
amplitude of optically bright quasars. Using the empirical 
correlation between black-hole mass and circular velocity of 
the host galaxy (Ferrarese 2002), in Section |1| we present a 
new derivation of the function N{M) based on the observed 
quasar luminosity function. We then show that this model 
is in agreement with the clustering measurements presented 
in Section |21 In Section |^ we present a Bayesian study of 
the halo occupation number which combines the results from 
Sectionsl^and^l This allows us to fully constrain all the pa- 
rameters of the halo model. Estimates of the quasar lifetime 
are presented in Section |S1 while the evolution of the halo 
occupation number over the redshift range 0.8 < 2 < 2.1 is 
addressed in Section Eventually, we discuss our results in 
Section |H| and conclude in Section |^ 

Throughout this work, we assume that the mass den- 
sity parameter flo = 0.3 (with a baryonic contribution 
fib ~ 0.049), the vacuum energy density parameter Oa — 0.7 
and the present-day value of the Hubble constant Hq — 
100/ikms"^Mpc"^ with h = 0.7. We also adopt a cold 
dark matter power spectrum with primordial spectral index 
n = 1 and ag = 0.8 (with erg the rms linear density fiuctu- 
ation within a sphere with a radius of 8h~^ Mpc). This is 
consistent with the most recent joint-analyses of tempera- 
ture anisotropies in the cosmic microwave background and 
galaxy clustering (see e.g. Tegmark et al. 2004 and refer- 
ences therein) and with the observed quasar power spectrum 
(Outram et al. 2003). 



2 THE DATA 

In this section, we present the main properties of our dataset 
and compute its clustering properties as a function of red- 
shift. 



2.1 Quasar selection and sample definitions 

The 2dF QSO Redshift Survey (2QZ) is a homogeneous 
database containing the spectra of 44,576 stellar-like objects 
with 18.25 < 6j < 20.85 (Groom et al. 2004). Selection of 
the quasar candidates is based on broadband colours (ufejr) 
from APM measurements of UK Schmidt Telescope (UKST) 
photographic plates. Spectroscopic observations with the 
2dF instrument (a multi-fibre spectrograph) at the Anglo- 
Australian Telescope have been used to determine the in- 
trinsic nature of the sources. The full survey includes 23,338 
quasars (the vast majority of which are endowed with a high- 
quality identification and/or redshift) which span a wide 
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Table 1. Main properties of our datasets. The superscripts min, max and med respectively denote the minimum, the maximum and the 
median value of a variable. 
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redshift range (0.3 ^ z < 2.9) and are spread over 721.6 
deg^ on the sky (see Groom et al. 2004 for further details). 

The 2QZ is affected by incompleteness in a number of 
different ways (for a detailed discussion see Groom et al. 
2004 ). In order to minimise systematic effects, we restrict 
our analysis to a sub-sample of sources defined by a min- 
imum spectroscopic sector completeness of 70 per cent. ^ 
Moreover, we only consider regions of the 2QZ catalogue for 
which the photometric completeness is greater than 90 per 
cent; this corresponds to a redshift range 0.5 < z < 2.1. 
Finally, we impose a cut in absolute magnitude, so that we 
only consider quasars brighter than Mtj — Slog^Q h — —21.7, 
which, assuming h = 0.7, is equivalent to Alb, = —22.5. Such 
an absolute magnitude cut ensures the exclusion of quasars 
where the contribution from the host galaxy may have led 
to a mis-identification of the source. 

In order to detect possible evolutionary effects, we want 
to subdivide our sample into three redshift bins. In partic- 
ular, we require that: (i) a similar number of quasars lies in 
each redshift bin, and (ii) each sub-sample covers a not too 
different interval of cosmic time. For this reason, we revise 
our initial sample selection by imposing an additional red- 
shift cut, so to keep only objects within 0.8 < z < 2.1. 
In fact, the time interval covered by the redshift range 
0.5 < z < 0.8 (1.78 Gyr) corresponds to one third of the 
total time elapsed between z = 2.1 and z = 0.5 (5.35 Gyr), 
whereas the number of quasars with 0.5 < z < 0.8 repre- 
sents less than 10 per cent of the selected quasar sample. 
By restricting the analysis to 0.8 < z < 2.1, we can greatly 
improve on both the previously mentioned conditions. More- 
over, we obtain a sample for which the mean number density 
of quasars is very weakly varying with redshift (as it can be 
seen in Fig. 1 of Outram et al. 2003), since through this cut 
we remove the largest mean number density variations as a 
function of redshift. The drop in mean density is less than 
60 per cent between z = 0.8 and z = 2.1. 

With the above selection, we end up with nearly 14,000 
quasars (split in two regions - the north galactic pole strip 
(NGP) and the south galactic pole strip (SGP) - with re- 
spectively '^ 7800 and ~ 6100 quasars each), of which 75 per 
cent reside in regions with total completeness larger than 80 
per cent. In order to study the evolution of quasar cluster- 
ing, we divide this sample into three redshift slices. To satisfy 
our previously mentioned criteria, we end up choosing the 
following three intervals: 0.8 < z < 1.3, 1.3 < z < 1.7 and 
1.7 < z < 2.1, containing each between ~ 4300 and ~ 4900 

^ A sector is as a unique area on the sky defined by the intersec- 
tion of a number of circular 2dF fields. 



quasars (see Table 0. We note that the time covered be- 
tween z — 0.8 and z = 1.3 (1.91 Gyr) is nearly twice the 
time covered between z — 1.3 and z = 1.7 (0.97 Gyr); this is 
however unavoidable if we want to keep similar numbers of 
quasars in each redshift interval. As the sample is magnitude 
limited, each redshift interval will correspond to quasars of 
different intrinsic luminosities, a point we address further in 
Section Hm 



2.2 Number densities of quasars 

Groom et al. (2004) provide an analytical fit for the 6j 
quasar luminosity function, in the case of sources brighter 



than Mh 



51ogio/i > -21.7 and for QA < z < 2.1. They 



model the optical luminosity function as a double power law 
in luminosity which, as a function of magnitude (number of 
quasars per unit magnitude per unit volume), becomes 



$(Af6,,,z) = 
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0.4,3i(M(, -Af* ) 



+ 10 
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where the evolution is encoded in the redshift dependence 
of the characteristic magnitude M^. = M^. {z) — M^. (0) — 
1.08fcr(z) with t{z) the fractional look-back time (in units 
of the present age of the Universe) at redshift z and k a 
constant. A table with the best-fitting values for the pa- 
rameters /9i, 132, Mi*j (0) and k is provided by Groom et al. 
(2004), together with their statistical uncertainties. 

Equation ^ can be used to compute the selection func- 
tion of the 2QZ between Zmin < z < z^^^: 
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where M^ (z) and A/f^j [z) denote, respectively, the bright- 
est and faintest absolute magnitudes which are detectable 
at redshift z. These are obtained by using the K-correction 
from the composite quasar spectrum by Brotherton et al. 
(2001), and by assuming fixed apparent magnitude lim- 

20.85. The integration 
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limits in the denominator of equation Q are M" 
min(,„,„,,„,,)M6';j(z) and M^^;"" = max(,„^,„,,.„^^)M6fj(z). 
The comoving volume effectively surveyed is then given by 
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(3) 
where fis denotes the solid angle covered by the survey 
and |dT^/d«df2| is the Jacobian determinant of the transfor- 
mation between comoving and redshift-space coordinates. 
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which gives the comoving volume element per unit redshift 
and solid angle. We can then estimate the mean number 
density of quasars by writing 



*QSO 



nQSoiZmin, Zmax) — /^ 



,_ V'cff l^^min; ^max J 



(4) 



where Nqso and ^^ Wi are, respectively, the total number 
of observed quasars in the range Zmin < z < Zmax and its 
completeness weighted counterpart. Results obtained by ap- 
plying equation @ are listed in Table where we summa- 
rize the main properties of our samples. Two types of errors 
are quoted for nqso- Those listed first are determined by 
independently varying the four parameters which define the 
optical quasar luminosity function (i.e. /3i, P2, M^ (0) and 
fc) within their 1 a uncertainties as reported by Croom et 
al. (2004). ^ On top of this error, we quote a ~ 10 per cent 
uncertainty on nqso due to the large-scale distribution of 
quasars: as in Outram et al. (2003), we typically find a '^ 10 
per cent difference in the number counts between the SGP 
and NGP. As this difference appears to be nearly constant 
as function of redshift, we quote for all subsamples the same 
typical uncertainty due to large scale structure. We note 
that our determination of nqsq is independent of the nor- 
malisation constant $* appearing in equation Q. 

From the luminosity function we also compute the 
effective redshift of the different samples as Zcft = 
EtT° »» 2./EtT° ^^ (see TableEJ. 

2.3 Quasar projected correlation function 

The simplest statistic which can be used to quantify cluster- 
ing in the observed quasar distribution is the 2-point corre- 
lation function in redshift space, ^°^^{r±,TT). This measures 
the excess probability over random to find a quasar pair sep- 
arated by n along the line of sight and by r± in the plane 
of the sky. These separations are generally derived from the 
redshift and the angular position of each source, so that the 
inferred tt includes a contribution from peculiar velocities. In 
consequence, the reconstructed clustering pattern in comov- 
ing space comes out to be a distorted representation of the 
real one and ^°'^''(rx, tt) is found to be anisotropic. To avoid 
this effect, and determine the quasar clustering amplitude in 
real space, one can then use the "projected correlation func- 



tion" which is obtained by integrating ^° 
direction: 



{r±, tt) in the tt 



' = -— / ^ (rx,7r)d7r. 



r± 



(5) 



and it is therefore unaffected by redshift-space distortions. 
In this section, we measure this quantity for our quasar sam- 
ples. 

We start by building a catalogue of unclustered points 
which has the same angular and radial selection function as 
the data. The angular selection is trivially given by the dif- 
ferent completeness masks (see Croom et al. 2004 for further 



details), and we modulate the number of random points laid 
down as a function of spectroscopic and photometric com- 
pleteness. The radial selection function is instead obtained 
by heavily smoothing the observed quasar redshift distribu- 
tion, M{z), or even the observed quasar comoving distance 
distribution, J\f{r). Both uniform and Gaussian smoothings, 
with characteristic smoothing length of several hundreds of 
h~^ Mpc, have been used. The quoted results are insensitive 
to the precise details of the modelling of the radial selection 
function. This is partially due to the fact that the volume 
covered by the quasar sample is very large and that there 
are not many clear groups or clusters of quasars: the quasar 
redshift distribution is rather smooth when compared to a 
standard galaxy distribution (e.g. Fig. 1 of Outram et al. 
2003 vs Fig. 13 of Norberg et al. 2002). 

The quasar correlation function is then estimated by 
comparing the probability distribution of quasar and ran- 
dom pairs on a two-dimensional grid of separations {r±,n). 
We use both the Landy-Szalay estimator (Landy & Szalay 
1993) and the Hamilton estimator (Hamilton 1993): 

,obs _ DD- 2DR + RR ^,, _ DDjRR 

^^^ - RR ' ^f^ - (DRY " ^ ^^^ 

where DD, DR and RR are the suitably normalised numbers 
of weighted data-data, data-random and random-random 
pairs in each bin. ^ As expected, the two estimators give 
comparable answers within the errors. For this reason, in 
what follows we only present results obtained with the 
Hamilton estimator. 

With the current quasar sample, we find that a reliable 
measure of ^°^''{r_i_,-K) is only achievable on scales tt < 50 
h~^ Mpc. In fact, the number of quasar pairs with small 
transverse separations and large line-of-sight separation is 
very small. This is partially due to the sparseness of the 
samples considered but also to the rareness of large struc- 
tures of quasars. Eventually, we compute the projected cor- 
relation function using equation Q. In order to avoid the 
measured signal to be dominated by noise, we limit the in- 
tegration to an upper limit, vTmax. This limiting value needs 
to be sufficiently large in order to give a reliable and mean- 
ingful measurement of E°''°(rx)/rx on the scales of interest 
(i.e. rx ^ 20 h~^ Mpc) but also small enough to be less 
sensitive to noise. Using the redshift distortion models by 
Hoyle et al. (2002), we find that Tr^ax = 45 h~^ Mpc ful- 
fils both requirements. All the data presented in this paper 
are obtained using this value. In any case, we verify that 
our results are not sensitive to the precise value adopted for 



2.4 Error estimates for clustering measurements 

It is common practice to estimate errors on the clustering 
measurements from mock surveys based on galaxy forma- 
tion models. However, for quasars, such catalogues are not 
publicly available. We therefore opt for a bootstrap resam- 
pling technique. For each redshift-interval we divide both 



^ If the errors for the four parameters are statistically indepen- 
dent, the quoted values for Anqgq approximately give 1 a un- 
certainties, whereas, if the parameters are correlated (which they 
most certainly are), the quoted errors for nqgo correspond to a 
higher confidence interval. 



^ Note that, in this case, there is no need to use the standard J3 
(minimum variance) weighting scheme (Efstathiou 1988) since the 
mean density of quasars, nqgq, is so low that 1 + 4 tt J3 nqsq ~ 1 
for any reasonable quasar clustering amplitude. 
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Figure 1. Left panel: The projected correlation function for different samples of quasars from the 2QZ (data with errorbars). For each 
redshift sub-sample, the best-fitting power law is represented with a long-dashed line. For reference, the best-fitting power law for the full 
sample (0.8 < z < 2.1) is shown with a short-dashed line. The continuous lines represent the best-fitting constant-bias models discussed 
in Section l2.YI These functions practically coincide with the best-fitting halo-occupation models presented in Sections 13.41 and ISl Data in 
the shaded regions are derived from less than 20 quasar pairs and are not considered for model fitting. Right panel: Contour levels for 
the likelihood function obtained by fitting the data with a power law model. The best-fitting models are marked with a dot and the lines 
correspond to 4 different values of Ax^ = X^ ~ Xmin' ^^ particular, for Gaussian errors, the inner contours (Ax^ = 1 and 2.3) mark the 
68.3 per cent confidence levels for one and two parameters, respectively. Similarly, the outer contours (Ax^ = 4 and 6.17) correspond to 
the 95.4 per cent confidence levels for one and two parameters. 



the NGP and the SGP samples into eight equal-volume re- 
gions, and we measure the clustering properties of each sub- 
sample. For ease of calculation, the division is just based on 
redshift and is such that the depth of each region is larger 
than the adopted value for Ttmax- Since the number density 
is roughly constant as function of redshift, each of these re- 
gions approximately contains the same number of quasars. 
We then build 1000 bootstrap-samples, each of them com- 
posed by sixteen sub-samples (eight for each strip) randomly 
drawn (allowing repetitions) from the set described above. 
We measure the projected correlation function for each arti- 
ficial sample by appropriately averaging over the number of 
quasar and random pairs of the sub-samples (and not over 
individual quasar clustering measurements) . For each r± , we 
identify the rms variation of H over the bootstrap-samples 
with the 1 a error for the projected correlation function. 



Our method for estimating errors relies on the fact that 
our dataset is statistically representative of the quasar dis- 
tribution in the Universe. However, this cannot be true for 
bins of spatial separations which contain just a few quasar 
pairs. Therefore, in what follows, we ignore clustering re- 
sults obtained by less than 20 quasar pairs. Depending on 
the sample, this corresponds to r^ < 1 — 2 h~^ Mpc. Note 
that, on such scales, corrections for close pair incomplete- 
ness due to fibre collisions (Groom et al. 2004) should be 
also taken into account. 



2.5 Results 

The projected correlation function we obtained for the dif- 
ferent redshift bins (and for the total sample) is presented in 
the left panel of Fig. A clear evolutionary trend emerges 
from the data: the clustering amplitude for the high-redshift 
sample is nearly a factor of 2 (3) higher than for the total 
sample (low-redshift one). 

As a zero-th order approximation, we fit our results 
with a power-law functional form. This phenomenological 
description has been commonly used in the literature and 
allows us to compare our results with previous studies. More 
detailed modelling is presented in the next sections. Here, we 
assume that the quasar 2-point correlation function scales 
as 



ar) 



ro 



(7) 



where r denotes the comoving separation between quasar 
pairs. The corresponding projected correlation function is 
obtained through the simple integral relation. 



:(rx) = 2 r 

J r I 
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which, in the power-law case, reduces to 

r(i/2)r[(7-i)/2] 



= r_L 



dr 



ro 



r(7/2) 



(8) 



(9) 



with r(a;) the Euler's Gamma function. We use a minimum 
least-squares method (corresponding to a maximum likeli- 
hood method in the case of Gaussian errors) to determine 
the values of ro and 7 that best describe the clustering data. 
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Table 2. Best-fitting power law and constant-bias models for the four quasar samples. The goodness of each fit is measured by the 
quantity Xmin/'^'-'f which gives the minimum value assumed by the chi-square statistic divided by the number of degrees of freedom. The 
symbol M[, denotes the halo mass which matches the observed clustering amplitude. 



''O 



(h-i Mpc) 



dof 



(h'^Mpc) 



dof 



^0 



(h-^ Mpc) 



dof 



dof 



A4 



0.8 


1.3 


5.4; 


1.3 


1.7 


4.3; 


1.7 


2.1 


7.6 



-1-0.9 
-1.3 
-1-1.8 
-2.0 
-1-1.2 
-2.1 



2.02 
1.49 
1.79 



+0.36 

0.33 

+0.32 

-0.35 

+0.25 
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0.7 
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2.85/3 
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0.7 

+0.4 
0.5 

+0.8 
'0.8 
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0.18 
0.19 
+0.32 
-0.35 
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0.84/3 13.26 



i2mtUl 



+0.11 
0.12 

+0.11 



0.8 



2.1 



D+0.9 
'-1.5 



1.53 



+0.19 
0.21 



0.13/2 



a+0.6 
'-0.6 



0.13/3 



5.4 



+0.5 



2.54/3 2.42 



+0.20 
0.21 



0.57/3 12.9llole 



A principal component analysis (see e.g. Porciani & Gi- 
avalisco 2002) is used here to deal with correlated errorbars. 
The principal components of the errors have been computed 
by diagonalizing the covaxiance matrix obtained by resam- 
pling the data with the bootstrap method described in the 
previous section. The objective function (the usual y^ statis- 
tic) has been obtained by considering the first four principal 
components which, for each redshift interval, account for 
more than 85 per cent of the variance. 

The best-fitting values for 7 and ro are given in Table 
|5|and the corresponding projected correlation functions are 
overplotted to the data in the left panel of Fig. Q Contour 
plots of the x^ functions are shown in the right panel of Fig. 
Note that the correlation length, ro, and the slope of the 
correlation function, 7, are strongly covariant: in order to 
accurately describe the data, larger values of ro need to be 
associated with steeper slopes. The best-fitting slope for the 



whole quasar sample, 7 = 1.53 



+0.19 



is in very good agree- 



ment with the redshift-space analysis by Croom et al. (2001) 
who found 7 = 1.56 



+0.10 



at a mean redshift of 



1.49. 



This is not surprising, since we only consider scales that are 
in the quasi-linear and linear regime of gravitational insta- 
bility where the correlation function in real and redshift- 
space are proportional to each other (Kaiser 1987). Accord- 
ingly, the comoving correlation length we find in real space, 
ro = 4.8l5'5 ^~^ Mpc, is, as expected, slightly smaller than 
its redshift-space counterpart, 5.691o;5q h~^ Mpc. '' 

As previously discussed, visual inspection of Fig. 0sug- 
gests that the quasar clustering amplitude at 1.7 < z < 2.1 
is nearly a factor of 2 higher than that obtained for the whole 
sample. Two sigma evidence for an increase in the cluster- 
ing amplitude of optically selected quasar between z ~ 1 and 
z ~ 1.8 has been presented by La Franca, Andreani & Cris- 
tiani (1998). However, given their sample size (a few hundred 
quasars) it is not clear whether the detected evolution is real 
or it is spuriously created by cosmic variance effects (see e.g. 
the discussion in Croom et al. 2001 who, using a preliminary 
data release of the 2QZ, found a that the redshift-space clus- 
tering amplitude at z = 2.7 is a factor 1.4 higher than at 
z = 0.7 which is marginally significant). It is therefore in- 
teresting to try to quantify the evolution of the clustering 

* Croom et al. (2001) assume statistically independent Poisson 
errorbars for the correlation function at different separations. This 
explains the factor of 2 between their and our estimate of the 
uncertainty for ro and 7. 



amplitude in our large quasar sample. In order to facilitate 
the comparison among the different subsamples (and with 
previous studies), we report in Table |21 the 68.3 per cent 
confidence intervals for ro obtained by assuming 7 = 1.53 



(r'o 



(7=1.53) 



) and 7 = 1.80 (r^' 



(7=1.80) X 



When fixing the slope. 



we note a steady increase of the quasar correlation length 
with z. Within this approximation, clustering evolution with 
redshift is detected at the ~ 3.6 it confidence level. However, 
since we are dealing with a flux-limited sample (where the 
highest-redshift objects have on average the highest intrinsic 
luminosities) , it is not possible to say whether this evolution 
of ro corresponds to a real change in the quasar population 
or to clustering segregation with luminosity. By analysing 
quasar clustering in the range 0.3 < z < 2.9 as a function 
of apparent luminosity in the preliminary data release cat- 
alogue of the 2QZ, Croom et al. (2002) found weak (~ 1 a) 
evidence for the brightest third of quasars on the sky to 
be more clustered than the full data set. Even though the 
different selection criteria prevent a direct comparison, we 
find a statistically more significant change in the cluster- 
ing strength among our sub-samples than that reported in 
Croom et al. (2002) . We defer a detailed analysis of the lu- 
minosity dependence of quasar clustering to future work. 



2.6 QSO vs galaxy clustering 

How do our results on the spatial distribution of quasars 
compare with galaxy clustering at similar redshifts? Un- 
til very recently, only rather small samples of high-redshift 
galaxies were available and any attempt to determine their 
clustering properties was most probably hampered by cos- 
mic variance (e.g. Le Fevre et al. 1996; Carlberg et al. 1997; 
Arnouts et al. 1999; Magliocchetti & Maddox 1999). The 
advent of color-selected surveys has allowed the detection 
of rich and homogeneous samples of high-redshift galaxies 
over relatively large volumes. We want to compare the re- 
sults obtained from the largest samples presently available 
with those obtained from our sample of quasars. A num- 
ber of studies have shown that Lyman-break galaxies at 
2: ~ 3 are strongly clustered (e.g. Porciani & Giavalisco 
2002 and references therein). Both their correlation length, 
ro ~ 4 /i~^ Mpc, and the slope of ^, 7 ~ 1.5, (Porciani & 
Giavalisco 2002; Adelberger et al. 2003) are remarkably sim- 
ilar to the values obtained from our quasar sample. On the 
other hand, star-forming galaxies at 2: ~ 1 (detected by ex- 
ploiting the Balmer break in their spectra) are found to be 
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Figure 2. Quasar-to-mass bias function derived by applying equation 1101 to quasar samples within different redshift ranges (points 
with errorbars). The continuous line shows the best fitting constant value. Dashed lines indicate the values of the bias for which Ax^ = 1 
(short dashed) and Ax^ = 4 (long dashed). 



slightly less clustered: ro ~ 3 h~^ Mpc with 'y ^ 1.7 (Adel- 
berger 2000). Similarly, the galaxy-galaxy correlation func- 
tion from the DEEP2 Galaxy Redshift Survey at z^s = 1.14 
is well described by a power law with 7 = 1.66 ± 0.12 and 
ro = 3.1 ± 0.7 /i"^Mpc (Coil et al. 2004). Evidence that 
early-type galaxies are more clustered than late-type ones 
(Firth et al. 2002; Coil et al. 2004) might help reconciling 
these results with those inferred from the quasar popula- 
tion. Extremely red galaxies at z ~ 1 have been found to 
be exceptionally strongly clustered. Assuming 7 = 1.8 (as 
inferred from their angular clustering) one finds ro ~ 12 
h-^ Mpc (Daddi et al. 2001; Firth et al. 2002; Roche, Dun- 
lop & Almaini 2003). 



2.7 QSO vs dark-matter clustering 

To study how the spatial distribution of quasars relate to the 
underlying matter distribution, we introduce the quasar-to- 
mass bias function: 



^('■X,2cff) = 



H(rx,2:„ 



< Z < Z-a 







Hm(rx,Zcff) 



1/2 



(10) 



where Hm is the projected correlation function of the mass 
distribution in the assumed cosmology computed as in Pea- 
cock k, Dodds (1996). Fig.|21shows our results for the differ- 
ent redshift bins. Errorbars at different spatial separations 
are not statistically independent. As previously described, 
we fit the data with a constant function by using the prin- 
cipal components of the bootstrap errors (which shows that 
datapoints at rx > 10 /i~^ Mpc are strongly correlated). 
Results and the corresponding 1 a uncertainties are listed in 
Table|21 We find that h steadily increases with z. This statis- 
tically significant trend is mostly due to the rapid evolution 
of the mass autocorrelation function. 

It is interesting to determine the halo mass, Afb, which 
corresponds to the observed quasar clustering amplitude (i.e. 
as if all quasars would reside in haloes with a fixed mass). 
We find that, for all the subsamples, Mb is of the order of 
10^^ Mq and that it slightly increases with z (see Table |5J. 



3 THE HALO OCCUPATION NUMBER OF 
QUASARS 

In this section, we present the halo model for the spatial 
distribution of quasars. After introducing the basic nota- 
tion, we describe the main features which characterize our 
model. ^ We then use the number density and the projected 
correlation function determined in Section |21 to constrain 
the free parameters of the halo model. This allows us to 
determine the way quasars populate dark-matter haloes of 
different masses. 



3.1 The halo model 

It is generally believed that quasars are powered by mass 
accretion onto supermassive black holes lying at the centre 
of galaxies. Cold dark matter (CDM) models for structure 
formation predict that galaxies form within extended dark- 
matter haloes. The number density and clustering properties 
of these haloes can be easily computed, at any redshift, by 
means of a set of analytical tools which have been tested 
and calibrated against numerical simulations (e.g. Mo & 
White 1996; Sheth & Tormen 1999). In consequence, the 
problem of discussing the abundance and spatial distribu- 
tion of quasars can be reduced to studying how they pop- 
ulate their host haloes. The basic quantity here is the halo 
occupation distribution function, Pjv(Af), which gives the 
probability of finding A'^ quasars within a single halo as a 
function of the halo mass, M. Given the halo mass function 
n{M) (number of dark-matter haloes per unit mass and vol- 
ume), the mean value of the halo occupation distribution 
N[M) = {N){M) = J2n^ Pn{M) (which from now on wiU 
be referred to as the halo occupation number) completely 
determines the mean comoving number density of quasars: 



/ 



n{M) N{M) AM 



(11) 



^ Further details can be found in Magliocchetti & Porciani (2003) 
where we used a similar tool to study the clustering properties of 
galaxies with different spectral types in the 2dF Galaxy Redshift 
Survey (2dFGRS). 
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Analogous relations, involving higher-order moments of 
Pjv(M), can be used to derive the clustering properties of 
quasars in the halo model. For instance, the 2-point corre- 
lation function can be written as the sum of two terms 



ar)=e\r)+e\r) 



(12) 



The function ^^ , which accounts for pairs of quasars re- 
siding within the same halo, depends on the second facto- 
rial moment of the halo occupation distribution, E^ (Af) = 
{N(N — 1))(M) and on the spatial distribution of quasars 
within their host haloes, p(x; M) , ® 



rw^/ 



n(M)S2(M) 



dM / p{x-M)p(x + r-M)d^x . (13) 



On the other hand, the contribution to the correlation com- 
ing from quasars in different haloes, ^qsq, depends on N{M) 
and p(x; M) as follows 



e\r)=i!<MimMiidM, 



n(M2)N(M2) 



dMa 



X / p(xi;Afi)p(x2;A/2)eh(ri2;Mi,M2)d\idV2 ,(14) 

where Yi marks the position of the centre of each halo, 
ri2 = r2 — ri is the separation between the haloes, Xi de- 
notes the quasar position with respect to each halo centre 
(hence r = ri2 -I- X2 — xi), and ^h(''i2; Afi, Af2) is the cross- 
correlation function of haloes of mass Afi and A^2, separated 
by ri2 . For separations which are larger than the virial radius 
of the typical quasar-host halo, the 2-halo term dominates 
the correlation function. In this regime, 5h(r; M\, M2) scales 
proportionally to the mass autocorrelation function, ^m(r), 
as Ch(»-; Mi,M2) ~ b{Mi) h(M2) ^m(r), with b{M) the linear- 
bias factor of haloes of mass M (Cole & Kaiser 1989; Mo & 
White 1996; Catelan et al. 1997; Porciani et al. 1998). As a 
consequence of this, the large-scale behaviour of the quasar 
correlation function also comes out to be ^(r) ~ ^cffCm(f), 
with 



/ b{M) N{M) n{M) AM 
/ N{M) n{M) AM 



(15) 



Note that all the different quantities introduced in this 
section depend on the redshift z, even though we have not 
made it explicit in the equations. 

In order to use the halo model to study quasar clus- 
tering, one has to specify a number of functions describing 
the statistical properties of the population of dark-matter 
haloes. In general, these have either been obtained analyti- 
cally and then calibrated against N-body simulations, or di- 
rectly extracted from numerical experiments. For the mass 
function and the linear bias factor of dark-matter haloes 
we adopt here the model by Sheth & Tormen (1999). We 
approximate the 2-point correlation function of dark-matter 
haloes with the function (see e.g. Porciani & Giavalisco 2002; 
Magliocchetti & Porciani 2003) 

,( M M\ JKAh)b{M2)^^{r) ifr>ri+r2 
^h(r;Mi,M2) = < (16) 

1—1 it r < ri + r2 



where the mass autocorrelation function, S,m{r), is computed 
using the method introduced by Peacock & Dodds (1996) 
which, for our purposes, is sufficiently accurate both in the 
linear and non-linear regimes. ^ For small separations, equa- 
tion II6II accounts for spatial exclusion between haloes (e.g. 
2 haloes cannot occupy the same volume). In Section 13.61 
where we discuss the small-scale clustering of quasars, we 
identify the Eulerian zone of exclusion of a given halo, ri, 
with its virial radius. At the same time, we assume that 
quasars trace the dark-matter distribution and adopt, for 
p{x;M) a Navarro, Frenk & White (1997; NFW) profile 
with concentration parameter obtained from equations (9) 
and (13) of Bullock et al. (2001). Note that, in all the other 
sections of this paper, we only consider the large scale distri- 
bution of quasars (r > 2 h~^ Mpc), where exclusion effects 
and density profiles do not affect the predictions of the halo 
model for ^. 



3.2 Clustering on the light-cone 

Equations II12L 11311 and 1141 1 describe the clustering prop- 
erties of a population of cosmic objects selected in a three- 
dimensional spatial section taken at constant cosmic time in 
the synchronous gauge. However, deep surveys like the 2QZ 
span a wide interval of lookback times and the equations we 
have introduced above do not apply in this case. 

A number of authors have discussed 2-point statistics of 
objects lying on the light-cone of the observer (e.g. Matar- 
rese et al. 1997, Yamamoto & Suto 1999, Moscardini et 
al. 2000 and references therein). These works have shown 
that the observed correlation function can be written as the 
weighted average: 



r»->V(2)^(r,z)d2 



J^.„.x y^;(^) j^ 



(17) 



where W{z) = J\f{zf{dV/dz)~'^ , with M{z) the actual red- 
shift distribution of the objects in the sample and dV/dz 
the Jacobian between comoving volume and redshift. Note 
that equation 11711 only holds for scales r over which: i) M is 
nearly constant and ii) ^ does not significantly evolve over 
the time r/[(l -I- z)c\ (where c denotes the speed of light 
in vacuum) . Within the range of separations covered by our 
dataset, both the conditions are satisfied for our quasar sam- 
ple. 

Combining equations (|HJ, 11211 and 11711 . we compute 
H°'^''(r-x) in our 4 redshift intervals for a large set of N{M) 
models and then compare the results with 3{r±,ZaH), the 
constant-time correlation function evaluated by using equa- 
tions JHJl and 11211 at the effective redshift of each sub- 
sample. In all cases we find extremely good agreement be- 
tween the two functions. Even for the widest redshift bin, 
0.8 < z < 2.1, we find a meiximum discrepancy of 2 per cent 
which is negligibly small with respect to the typical error as- 
sociated with the observed correlation function. Therefore, 
in what follows, we use H(rx, Zch) to compare the predictions 
of different models with the clustering data. This greatly 



^ This is normalised in such a way that Jn^" p(y; M)cPy = 1 
where Tvir is the virial radius which is assumed to mark the outer 
boundary of the halo. 



^ In principle, for separations where £,m.{r) ~ 1, non-linear terms 
should be added to equation 1161 (Mo & White 1996; Catelan 
et al. 1997; Porciani et al. 1998). However, for the haloes and 
redshifts of interest, these can be safely neglected. 
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simplifies (and speeds-up) the model fitting procedure. As 
an additional test, the comparison between H° "(rx) and 
S(rx,Zoff) is repeated for all the best-fitting models that 
are discussed in the forthcoming sections and no significant 
difference is found. 



3.3 The halo occupation number 

The final, key ingredient needed to describe the clustering 
properties of quasars is their halo occupation distribution 
function. In the most general case, Pjv(M) is entirely spec- 
ified by all its moments which, in principle, could be obser- 
vationally determined by studying quasar clustering at any 
order. Regrettably, as we have already shown in Section|53| 
quasars are so rare that their 2-point function is already 
very poorly determined so that it is not possible to accu- 
rately measure higher-order statistics. As in Magliocchetti 
& Porciani (2003), we overcome this problem by assuming 
a predefined functional form for the lowest-order moments 
of Pn{M). It is, in fact, convenient to describe the halo oc- 
cupation number, N{AI), and (if necessary) its associated 
scatter, E'^(M), in terms of a few parameters whose values 
will then be constrained by the data. 

We consider here the "censored" power-law model. 



iV(M)=iVo(- 



e(M - Mo) 



(18) 



(where Q{x) denotes the Heavyside probability distribution 
function) which has been widely used in the literature to de- 
scribe galaxy clustering (e.g. Magliocchetti & Porciani 2003 
and references therein). In this case, the halo occupation 
number vanishes for M < Mq and scales as a power law of 
the halo mass for M > A'lo- The parameter A'^o gives the 
mean number of objects contained in a halo of mass Mq. * 
Studies of the local galaxy population with both hydrody- 
namical simulations and semianalytic models for galaxy for- 
mation are consistent with equation 1181 (e.g. Sheth & Di- 
aferio 2001; Berlind & Weinberg 2002; Berlind et al. 2003). 
We assume that the same parameterization is adequate for 
quasars at z > 1. Given the observational evidence for a cor- 
relation between black-hole and host-halo masses (Ferrarese 
2002; see also the discussion in Section FASII . it is reasonable 
to expect the presence of a threshold mass for the host haloes 
of a quasar population with a given minimum luminosity. At 
the same time, a power-law scaling with an unspecified in- 
dex Q for M > Mo is general enough to explore a wide range 
of possibilities. 

It would be ideal to test equation (1181 against physical 
models for quasar activity. A number of authors recently de- 
veloped simplified schemes to include black-hole accretion in 
galaxy formation models (e.g. Kauffmann & Haehnelt 2002; 
Enoki, Nagashima & Gouda 2003; Menci et al. 2003; Di 
Matteo et al. 2003 and references therein). Regrettably, at 
variance with studies of normal galaxies, mock catalogues 
produced with these models are not publicly available. For 
this reason, to test the reliability of equation (1181 . we are 
forced to follow an indirect approach by associating quasar 



* Note that equation 1181 is more general than the commonly 
used N{M) = (M/Nh)" e(A/ - Mq) which, for a = 0, automat- 
ically implies N{M) = 1 at any M > A/q- 



activity with a particular subset of synthetic galaxies. In 
particular, we consider galaxies which, at z > 1, contain a 
substantial amount of cold gas in their nuclear region that, 
in principle, could be accreted onto a central supermassive 
black hole. Thus, in Appendix B, we use a publicly available 
semianalytic model for galaxy formation (Hatton et al. 2003) 
to study the halo occupation number of galaxies which, at 
2; ~ 1, are actively forming stars in their bulges. The rea- 
son for selecting this sample is threefold: i) Even though 
imaged quasar hosts are consistent with being massive ellip- 
ticals (e.g. Dunlop et al. 2003), there is some observational 
evidence that, at high redshift, quasars might be associated 
with active star-formation (e.g. Omont et al. 2001; Hutch- 
ings et al. 2002). ii) In the local Universe, powerful Type 
2 active galactic nuclei are found in bulges with either on- 
going star formation or young stellar populations (Kauff- 
mann et al. 2003). iii) The observed correlation between 
black-hole mass and bulge velocity dispersion suggests that 
quasar activity and bulge formation probably are physically 
associated phenomena. For instance, galaxy interactions and 
bar-induced infiows might funnel some gas into the nuclear 
region of a galaxy thereby triggering simultaneous star for- 
mation and quasar activity. 

Our results show that the halo occupation number of 
the simulated galaxies is well approximated by a power-law 
with a rather sharp cutoff at low masses. This provides ad- 
ditional motivation to use equation 1181 in our analysis. 



3.4 Constraints from large-scale clustering 

Assuming equation 1181 , we apply a least-squares method to 
determine the values of Mq and a which best describe the 
clustering data presented in Fig. As discussed in Section 
12.51 we use a principal component analysis of the errors to 
deal with the clustering data. We only consider the region of 
parameter space where a > and Mo > 10^ Mq . We impose 
this lower limit to a since we expect the halo occupation 
number to be a non-decreasing function of the halo mass. ^ 
On the other hand, we assume a lower limit for Mo because 
it is unreasonable to consider halo masses which are smaller 
than the minimum inferred mass of the black holes powering 
our quasars (cf. Section IA4II . In the range of separations 
covered by our dataset (r > 2 h~^ Mpc), the two-halo term 
dominates the amplitude of the quasar 2-point correlation 
function and there is no need to specify the form of the 
function S^(Af). From equation 1151 . we also note that the 
quasar correlation function on large scales does not depend 
on the overall normalisation of N{M) (e.g. the parameter 
A^o). 

Contour plots of the x^ function are shown in Fig. |5] 
Note that, independently of the redshift interval, the region 
of parameter space allowed by the data is rather large and 
does not provide tight constraints on the values of Mq and 
Q. This is because our data only fix the normalisation of the 



Note that solutions with q < are allowed by the data if 



Mo ~ 10^ 



Mq (the precise value slightly depending on the 



assumed redshift range). In this case, quasars are hosted by haloes 
lying in a narrow mass range which is centred around the values of 
Mj, given in Tablel^ For this reason, there is no need to rediscuss 
these solutions here. 
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Figure 3. Contour levels of the x^ function for the parameters a 
and Mq obtained by fitting the large-scale clustering of quasars in 
different redshift intervals with the predictions of the halo model 
given in equation lliji . The various panels contain contour plots 
for the x^ surface in the a — Mq plane. The contours correspond 



to Ax"' 



■X 



1 and 4 (respectively marking the 68.3 and 



95.4 per cent confidence levels for two fully degenerate Gaussian 
variables). Contours of the value of logj^Q A^o which matches the 
observed number density are plotted as a function of a and Mq 
(dashed lines with labels) for the different redshift bins consid- 
ered. The dotted lines mark the halo masses which correspond to 
the observed clustering amplitudes (see Table 1^ . 



correlation function (e.g. the bias parameter shown in Fig. 
^ and there is a whole one-dimensional family of (Qf,A/o) 
pairs which correspond to the same clustering amplitude. 

3.5 Constraints from the number density 

In the halo- model described by equation (I18II . the number 
density of quasars depends on all the free parameters a, 
A/o and A^o- In particular, A'o acts as an overall normalising 
factor so that, for given values of a and A/q, it can be ex- 
pressed in terms of the mean number density of quasars by 
combining equations IIH and 1181 . This gives an additional 
relationship among the parameters of the model as, for ev- 
ery (q, Afo) pair, there is always a value of A'o which exactly 
matches nqso- The dashed lines in Fig. |3] show the regions 
in the {a, Mo) plane where the observed density corresponds 
to a given value of log^Q A'o (indicated by the labels in the 
plot). The parameters A'o and Mq are strongly covariant: in 
order to get the right quasar abundance, one needs to lower 
the normalisation of the halo occupation number when A/q 
is reduced. The allowed range for A'o spans many orders of 
magnitude, reflecting the steep slope of the halo mass func- 
tion at the low- mass end. 

3.6 Constraints from small-scale clustering 

For separations smaller than the typical size of the host 
haloes, the galaxy 2-point correlation function is dominated 



by the contribution of pairs lying within a single halo. In 
this regime, ^(r) is fully described by equation (I13t which 
encodes information on the halo occupation distribution 
through its second factorial moment, E^(Af). This function 
can be conveniently expressed in terms of the halo occupa- 
tion number as follows: 



E^(Af) =r{M)N{Mf 



(19) 



For a Poisson distribution, r{M) — 1 independently of M. 
In this case, measuring the 2-point correlation function on 
small scales provides additional constraints on the halo oc- 
cupation number. 

However, in general, the halo occupation distribution 
function is not Poissonian and T{M) depends on the halo 
mass. In principle, this complicates the estimate of N{M) 
from analyses of small-scale clustering. In fact, a number of 
additional free-parameters might be required to describe the 
behaviour of r(Af). On the other hand, though, models of 
galaxy formation suggest that, independently of the galaxy 
sample considered, r(Af) is a very simple function which 
can be parameterized in terms of the same variables that 
are used to describe N{M). Consistent results have been 
obtained for low-redshift galaxies by using different semi- 
analytical models (e.g. Sheth & Diaferio 2001; Berlind & 
Weinberg 2002) and hydrodynamical simulations (Berlind 
et al. 2003). Similarly, in Appendix B, we use a publicly 
available semianalytic model to study the function r(Af) 
for star- forming galaxies at 2: ~ 1. In all cases, the scatter 
of Pm(M) is strongly sub-Poissonian for haloes which, on 
average, contain less than 1 object, and nearly Poissonian 
for larger haloes. This property plays a fundamental role in 
breaking the degeneracy among all the models for the halo 
occupation number which can otherwise accurately describe 
galaxy clustering on large scales (Magliocchetti & Porciani 
2003). 

We do not know whether the same conclusions apply 
to quasars. It is anyway interesting to understand what this 
would iniply. Let us assume that, also for quasars, F <g 1 
when A' < 1 while F ~ 1 when N >l. Within the allowed 
parameter range in Fig.|3 this implies that, at variance with 
galaxies, the 1-halo term never dominates the quasar corre- 
lation function even on scales which are much smaller than 
the typical halo size. This happens because the quasar N{M) 
is always much smaller than unity and its associated scatter 
is strongly sub-Poissonian. In other words, the distribution 
of (optically bright) quasars in a halo is binary (either there 
is one or there is none) and it is basically impossible to find 
2 quasars being hosted by the same halo. 

From the absence of quasar multiplets in a single halo 
it follows that, in order to use observational determinations 
of quasar clustering on small scales to break the degeneracy 
among models which predict the same clustering amplitude 
on large scales, one has to rely on the detection of halo 
exclusion effects (different haloes cannot overlap) . Note that 
this would be a direct "measure" of the spatial dimension of 
dark-matter haloes and therefore of their mass. 

The exact signature induced by spatial exclusion is hard 
to predict since dark matter haloes are expected to be tri- 
axial objects and the precise form of the quasar correlation 
function on small scales is also expected to depend on the 
position of each quasar within its host halo (see equation 
I13II . However, it is clear that the configuration which max- 
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imizes this effect is obtained when quasars sit at the centre 
of the haloes. In this case, the distribution of quasars is a 
perfect (sparse sampled) tracer of the underlying halo distri- 
bution and the 2-point correlation function, £,{r), is expected 
to reach the value -1 on scales smaller than the typical size of 
the host halo. Such exclusion effects will then correspond to 
a flattening of the projected function H on the same scales. 
However, these effects might be hard to detect due to the 
small number statistics of close quasar pairs (cf. Fig. 0. 

Similar arguments apply to any population of rare ob- 
jects. A possible detection (with relatively low statistical 
significance) of exclusion effects has been reported from the 
analysis of the clustering properties of Lyman-break galaxies 
at redshift ~3 (Porciani & Giavalisco 2002). 

Note, however, that the scatter of the halo occupation 
distribution might depend on the detailed physical processes 
which give rise to the quasar phenomenon and thus be very 
different from the V function which describes the galaxy dis- 
tribution. Therefore, the presence of quasar multiplets inside 
single haloes is not ruled out by the data. Measuring the 
quasar 2-point correlation function on separations smaller 
than 1 h~^ Mpc would give a definitive answer to this ques- 
tion. 



4 THE HALO OCCUPATION NUMBER FROM 
QUASAR LUMINOSITIES 

Recent studies of stellar and gas dynamics in local galaxies 
have revealed a wealth of information on the population of 
supermassive black holes. The observational evidence for a 
correlation between the mass of a black hole, Mbh, and the 
circular velocity, Vc, of its host galaxy (Ferrarese 2002; Baes 
et al. 2003) is one of the most intriguing results. In this 
section, we use this empirically determined relation (Mbh oc 
Wc'^) to derive the halo occupation number of quasars in the 
2QZ. This is obtained by first converting quasar luminosities 
into a distribution of black-hole masses (to which we apply 
the i\fbh — Vc correlation) and then linking, with minimal 
assumptions, the circular velocity of the host galaxies to the 
mass of their dark-matter haloes. 



4.1 The mass of quasar host haloes 

The bolometric luminosity of a quasar and the mass of the 
accreting black hole can be related as follows 



Mbh _ f ibol 

Mq ^ Tj 1.26 X 10^8 erg s- 



(20) 



where 77 denotes the ratio between the bolometric luminosity 
of the quasar and the Eddington luminosity. We use this re- 
lation to determine the function V{M\Mb-f) which gives the 
conditional probability distribution of the host-halo mass, 
M, for a quasar with given absolute magnitude M^j . For ease 
of reading, we just summarize our calculations here while a 
detailed presentation of the model is given in Appendix A. 
To compute P(M|Mbj), we first determine the condi- 
tional probability of Mbh given A/^j . This is obtained by 
combining the most recent bolometric corrections from the 
B band with an empirically determined distribution of Ed- 
dington ratios (McLure & Dunlop 2004). We then use the 
observed Mbh — Vc relation (Baes et al. 2003) to derive the 



conditional probability of Vc given Mb , . Eventually, we con- 
vert circular velocities into halo masses by assuming that 
the observed circular velocity and the virial velocity of the 
host halo are related by Dc = V' '^vir with ^/j a free parame- 
ter. Recent lensing studies (Seljak 2002) suggest that, in the 
mass range of interest, %j) — 1.4 ±0.2 (case A); alternatively, 
theoretical arguments based on the estimated low concen- 
tration of high-redshift haloes, suggest that 1/) ~ 1 (case B). 
These different choices bracket the range of plausible values 
for ^ (see the detailed discussion in Appendix A). 



4.2 The halo occupation number 

In this section, we test whether the conditional mass distri- 
bution P(M|M(,j) is consistent with the quasar clustering 
data we measured in Section r2.5l By integrating over the lu- 
minosity function, 7'(M|Afi,,) can be easily turned into the 
mass function of dark haloes which are quasar hosts, 

n^{M,z)^ ■ $(Af6,,z)P(M|M6,)dM6, . (21) 

The corresponding multiplicity function, M nq{M,z), is 
shown in the left panel of Fig. 2] for different values of z. 
Independently of redshift, the halo mass distribution per 
unit logarithmic interval of M peaks around 1O^^'^~^'^M0: 
this is the characteristic mass of quasar hosts, whose value 
is comparable with those listed in Table |5| By then divid- 
ing nq{M,z) by the halo mass function, we obtain a new 
estimate for the halo occupation number 



N{M, z) 



nc^(M,z) 
n{M,z) 



(22) 



Note that this, in principle, could give rise to a biased es- 
timate of N{M, z). In fact, in a CDM cosmology, each viri- 
alized halo contains a number of sub-haloes within its rvir, 
and at least some of these sub-haloes will be associated with 
galaxies which formed within their local overdensities and 
then failed into the larger halo. Since the rotational prop- 
erties of galaxies are expected to be related to local dark- 
matter overdensities, the values we derived for M most likely 
refer to sub-haloes. On the other hand, n{M) describes the 
meiss distribution of the parent haloes. We can account for 
this problem by introducing the conditional mass function 
of the sub-haloes of mass Ms which lie within a parent halo 
of mass Mp > Ms, n(Ms|Mp). The probability that an halo 
with mass M is a parent one is then given by 



Pp(Af) 



n{M) 



n{M) + /j^7 n(Mp) n(M|Mp) dA/p ' 



(23) 



where the integral at the denominator gives the mass func- 
tion for the sub-haloes. We use two different functional forms 
for n(Ms|Afp) which have been derived from high-resolution 
numerical simulations (Sheth & Jain 2003; Vale & Ostriker 
2004 and references therein). In both cases, we find that 
the sub-halo correction is negligibly small. In fact, at the 
redshifts spun by our quasar sample, the haloes we are con- 
sidering are rather massive and only a few per cent of them 
have been included into larger units. 

Results for the halo occupation number are presented in 
Fig. El In all cases, for M < 10^* Mq, N[M) is weU approxi- 
mated by a broken power law which qualitatively resembles 
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Figure 4. Left panel: the multiplicity function (differential num- 
ber density per log-unit of mass and per unit volume) of haloes 
which are hosts of 2QZ quasars (continuous line for case A, long- 
dashed line for case B) is compared with the corresponding dis- 
tribution for all the dark-matter haloes (short-dashed line) at 
different redshifts. Right panel: corresponding halo occupation 
numbers obtained by taking the ratio of the quasar-host mass 
function to the total halo mass function (continuous line for case 
A, long-dashed line for case B). For comparison, the case A so- 
lution at 2 = 1.06 is reproduced with a short-dashed line in the 
two top panels. 



equation 11811 . The low-mass tail has a typical slope ~ — 3.6 
for case A and ~ —4.2 for case B, independently of redshift. 
On the other hand, the high-mass tail gets progressively 
steeper when moving from low to high redshifts. For case A, 
we get a flat N{M) for the low-redshift sample, to be com- 
pared with a slope of ~ 0.4 for the dataset at intermediate 
redshifts and of ~ 0.7 for the highest-redshift quasars. The 
corresponding numbers for case B are ~ 0.4, ~ 0.9, ~ 1.1. 
For M > lO^^M©, in both cases the halo occupation num- 
ber starts growing exponentially. This happens because, in 
the high- mass tail, nq(M, z) does not drop as fast as the ex- 
ponential cutoff of the halo mass function (see the left panel 
in Fig. |1| and the shaded region in the right panel). Most 
likely this is a spurious effect due to the simple assumptions 
we use to derive nq{Al,z). This artifact, however, does not 
affect our conclusions since the fraction of quasars that are 
found to reside in such massive haloes is extremely small (at 
most, a few xlO"'' for case A and less than 2 per cent for 
case B). 

Note that estimates for N{AI) based on equation 12211 
are obtained without any information on the clustering prop- 
erties of quasars. It is therefore interesting to check whether 
they are in agreement with the determination of H°'^''(rx) 
we presented in Section 12.51 For this reason, we compute 
the effective bias associated with the different halo occu- 
pation numbers presented in Fig. |1| For the lowest-redshift 
sample, we get bes = 1.63 (2.02) for case A (B). At inter- 
mediate redshifts, we find &cff = 2.14 (2.73) for case A (B). 
while, for the highest-redshift interval, we derive 6eff = 2.58 



(3.36) for case A (B). These numbers have to be compared 
with the observational results presented in Table |5| For case 
A, we note that, even though the estimated bias parame- 
ter increases with redshift (as like as the data in Table I^J, 
its value is in general too low to accurately describe the 
observed clustering. In other words, case A tends to under- 
estimate the mean mass of quasar host haloes. Predictions 
for case B, instead, are of better quality. In this case, our 
results for the bias parameter are rather accurate for the in- 
termediate redshift sample, while they tend to overestimate 
(underestimate) the observational results at low (high) red- 
shifts. Anyway, the bias inferred from our models is always 
acceptable (with respect to the statistical errors associated 
with the determinations of the projected correlation func- 
tion). The maximum discrepancy appears at Zcb = 1.89 and 
corresponds to a statistical significance of 1.4 ct. These re- 
sults are in agreement with the recent analysis by Wyithe & 
Loeb (2004) who showed that quasar models with Mbh oc v^, 
tjj — 1 and 77 ~ 0.1 — 1 are able to reproduce the evolution 
of the correlation length measured in a preliminary data re- 
lease of the 2QZ (Groom et al. 2001, 2002). More accurate 
clustering measurements are then required to detect possi- 
ble changes in the Mbh — Vc correlation and to distinguish 
them from effects due to evolution of other parameters (for 
instance rj, tp or the bolometric correction). 



5 A BAYESIAN ANALYSIS OF THE HALO 
OCCUPATION NUMBER 

We have shown that, at variance with galaxy clustering, lack 
of information on (and from) the 2-point correlation function 
of quasars on small scales does not allow us to break the 
degeneracy among the best-fitting models presented in Fig. 
121 In this section, we adopt a Bayesian approach and use 
information on quasar luminosities to further constrain the 
parameters of the halo occupation number. 

Assuming equation 11811 . for each redshift interval, we 
translate the probability density function for the errors of 
nQso and H(rx) into a likelihood function for the model 
parameters and we write 

£tot(a,Mo,iVo) = /:ciu.t(a,Mo)-£dons(a,Mo,iVo) . (24) 

where £ciust accounts for the large-scale clustering analy- 
sis presented in Section [3.41 and £dons for the quasar num- 
ber density. ^° In other words, the number density con- 
straints weights as much as a single independent point in 
the clustering analysis. We then apply Bayes' theorem to our 
dataset. For simplicity, to express our lack of prior knowl- 
edge, we adopt constant (non-informative) prior distribu- 
tions for logj^f, A'^o and a. On the other hand, as a prior 
distribution for logjQ(Mo/M0), we use the probability dis- 
tribution 7^(1/0 1 Af^^^) that we derived in Appendix A and 
presented in Section [4. II This is the probability distribution 
of the halo masses which harbour the faintest quasars that 
can be detected in each 2QZ sample considered. This prior 
knowledge is based on the empirically determined correla- 
tion between black-hole masses and the circular velocity of 



^^ Assuming Gaussian errors, —2\nCi 
denotes the usual chi-square statistic. 
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Figure 5. Posterior probability distribution for the parameters of the halo model. The four frames correspond to different redshift ranges 
as indicated by the labels. Thick and thin lines respectively refer to Prior A and B. Top-right: Contours of the joint distribution of a and 
Mq (obtained by marginalizing the three-dimensional posterior probability over A^o). The most-probable point is marked with a small 
circle. To facilitate the comparison with Fig.|2l the continuous lines show the points where —0.5 In Pmax/f = 1 and 4 (which, in this case, 
do not have any special meaning). As in Fig.|21 in order to represent the covariance of the different parameters, the dashed lines show 
the loci in the (a, Mq) plane where a given value of logjQ A^o (indicated by the labels) perfectly matches the observed number density 
of quasars. Other panels: The probability densities for each single parameter (obtained by marginalizing the posterior distribution over 
the remaining two variables) are shown in the top-left (for log^gMo/MQ), bottom-left (for logj^Q A^o) a-^d bottom-right (for a) panels. 
In the top-left panels, the dashed lines show the assumed prior distributions for Mq. 



the host galaxies. As previously discussed, to test the robust- 
ness of our method with respect to underlying systematic 
uncertainties, we consider two diflerent prior distributions 
corresponding to %p = 1.4 ± 0.2 (Prior A) and ^ = 1 (Prior 
B). 

Contours of the posterior probability in the a — Mo 
plane and the probability distribution of the single param- 
eters (marginalized over the remaining ones) are shown in 
Fig. 13 The corresponding best-fitting values and credibility 
intervals for the different parameters are listed in Table |21 

Note that adopting our informative prior on A/q is 
enough to break the degeneracy among the parameters of 



the best-fitting models. In practice, both priors exclude the 
region Mq < IO'^^Mq where haloes are too small to har- 
bour bright quasars. This is sufficient to determine a non- 
degenerate solution for each redshift range. The main char- 
acteristics of these solutions can be summarized as follows. 
In general, the cutoff mass. Mo, has a very mild evolution 
with redshift. Using prior A, we get Mo ~ (1 - 3) x IO^^Mq, 
while, with prior B, we obtain Mq ~ (2-5) x 10"Mq. " On 



^^ Note that Mq is only mildly covariant with the high-mass 
slope, a, in the sense that slightly lower values for Mq are gener- 
ally associated with larger values of a. 
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Table 3. Best-fitting values for the parameters of the halo model for different redshift ranges (superscript bf). The last three columns 
list the central 68.3 per cent credibility intervals for each single parameter. These have been obtained by marginalizing the posterior 
probability distribution function over the remaining parameters. The quoted values correspond to the 15.85, 50 and 84.15 percentiles. 
In the last column, we list the central 90 per cent credibility intervals for the quasar lifetime. We assume that a halo of 2 X 10^'^ Mq 
harbours, on average, a single super-massive black hole so that the halo occupation number of bright quasars coincides with their duty 
cycle. 
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the other hand, in order to match the rapidly evolving bias 
parameter of the three quasar samples with a nearly invari- 
ant Mo, the high-mass slope, a, tend to become steeper and 
steeper with increasing z. This is in qualitative agreement 
with the results presented in Section r4.2l At Zcff = 1.06, 1.51 
and 1.89, we respectively find a ~ 0.5, 0.8 and 1.4 for prior 

A, and a '^ 0.4, 0.6 and 1.1 for prior B. It is important to 
stress, however, that the parameter a is typically poorly de- 
termined. Strictly speaking, the data just set an upper limit 
for it. The allowed range for the normalisation parameter 
No varies systematically with the assumed prior. In brief. 
No spans a broader range (approximately from 3 x 10"'' to 
2x 10~2) when prior A is used. On the other hand, with prior 

B, the probability distribution for A^o is more peaked and 
ranges from 3 x 10""' to 3 x 10"^. We note that, for both 
priors, A^o is less tightly determined for our high-redshift 
sub-sample. 



6 THE QUASAR LIFETIME 

The number of optically bright quasars per halo can be 
used to estimate the duty cycle of quasar activity and, thus, 
the quasar lifetime (Haiman & Hui 2001; Martini & Wein- 
berg 2001). In brief, let us assume, for simplicity, that each 
dark-matter halo contains one supermassive black-hole. In 
this case, the fraction of active quasars per halo coincides 
with the quasar duty cycle. Assuming that quasar activity is 
randomly triggered (for instance by tidal interactions with 
neighbours) during the halo lifetime, in, the duty cycle can 
then be expressed as fq/in where fq denotes a characteris- 
tic timescale over which the quasar is visible in the optical 
band. Both a single optically bright phase and a series of 
shorter bursts correspond to the same fq. 



6.1 Estimating the quasar lifetime 

In this section, we use the posterior probability distribution 
presented in Fig. |3 to determine the characteristic quasar 



lifetime. For each halo mass, we first compute the probability 
density function of the halo occupation number: 



v{N)= / 5[N - N{M\e)]P{e\D)<r 



(25) 



where 5(x) denotes the Dirac-delta distribution, 6 = 
{a, Mo, No), P{9\'D) is the posterior probability and 
N{M;8) is given in equation ([TSl . The 5, 50 and 95 per- 
centiles of this distribution are shown in Fig.|S|as a function 
of the halo mass. Note that results are nearly independent 
of the considered prior. For haloes with M ~ IO'^^Mq, the 
occupation number is tightly constrained by the data. We 
thus estimate the characteristic quasar lifetime by assum- 
ing that, for these haloes, A'^ = tq/tu- Following Martini & 
Weinberg (2001), the halo lifetime is defined as the median 
time interval during which a halo of mass M at redshift z 
is incorporated into a halo of mass 2M. This quantity is 
computed using equation (2.22) of Lacey & Cole (1993). 

Results for tq are listed in Table |21 We find that 
the estimated quasar duty cycle increases with z (and/or 
with quasar luminosity). For our sample at ZcS ~ 1.06, 
we find that only ~ 0.6 per cent of the host-haloes with 
M = 2x IO'^^Mq contain a bright quasar, which corresponds 
to iq ~ 2 X 10^ yr. This coincides with the e-folding time of 
a black hole which accretes mass with a radiative efficiency 
e ^ 0.1 and shines at a fraction rj ~ 0.5 of its Eddington lu- 
minosity (Salpeter 1964). On the other hand, at Zcff — 1.89, 
the fraction of active black-holes increases to ~ 5 per cent 
and iq ~ 7 X 10^ yr. In all cases, the estimated lifetime lies 
between 10^ and 10** yr. 

Even though the determination of N{M) becomes more 
uncertain for M 3> W^^Mq, our results suggest that the oc- 
cupation number of quasars tends to increase with the halo 
mass. This, however, does not imply that fq augments as 
well. In fact, our estimates for the quasar lifetime are de- 
generate with the occupation number of supermassive black 
holes which, most likely, increases with the halo mass. 

Note that, given the simplicity of the model, our results 
are only indicative. The quoted quasar lifetimes should be 



Cosmic Evolution of Quasar Clustering 15 



revised upwards if: (i) a non-negligible number of haloes do 
not harbour any supemiassive black-hole; (ii) optical radia- 
tion from quasars turns out to be significantly beamed; (iii) 
in the presence of an important fraction of obscured sources. 
On the other hand, iq is smaller than what reported here if 
more than one supermassive black-hole is hosted, on aver- 
age, by each halo. 

A number of observations hint towards a one to one 
correspondence between supermassive black holes and host 
haloes. High-resolution optical imaging with the Hubble 
Space Telescope (HST) shows that bright quasars [My < 
—23) at z < 0.5 are only harboured by exceptionally lumi- 
nous galaxies with L > Ly (Bahcall et al. 1997; Hamilton, 
Casertano & Turnshek 2002) . These galaxies turn out to be a 
mixture of difi'erent morphological types, ranging from nor- 
mal ellipticals and spirals to complex systems of gravitation- 
ally interacting components (Bahcall et al. 1997). However, 
a number of observational results suggest that spheroidal 
hosts become more prevalent with increasing nuclear lumi- 
nosities: quasars with Mv < —23.5 are virtually all har- 
boured by luminous elliptical galaxies (Dunlop et al. 2003) . 
Similarly, observed surface-brightness profiles suggest that 
bright quasars at z ~ 1 — 2 are hosted by massive ellipticals 
undergoing passive evolution (Kukula et al. 2001; Falomo et 
al. 2004). ^^ Taking this for granted, we can show that our 
assumption of one supermassive black hole per halo (and 
thus our inferred quasar lifetime) is rather realistic. The ar- 
gument proceeds as follows, (i) Massive elliptical galaxies are 
made of old stellar populations which formed at 2; > 2 and 
passively evolved thereafter, (ii) In the assumed cosmology, 
the massive haloes which harbour these galaxies can only in- 
crease their mass by a factor of a few from z = 1 — 2 to the 
present epoch. This mainly happens via accretion of smaller 
objects, (iii) From clustering studies in the local Universe, 
we derive that haloes with M ~ 1Q^^~^'^ Mq harbour on av- 
erage ~ 1 — 2 early-type galaxies with effective luminosity 
Lcff = 1.3 L* (Magliocchetti & Porciani 2003). Points (i), (u) 
and (iii) imply that the mean number of early-type galaxies 
per halo was of order unity even at 2 ~ 1 — 2. Thus, our 
thesis follows from the assumption that each galaxy hosts a 
single supermassive black hole. 
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Figure 6. Halo occupation number obtained from the posterior 
probability distribution in Fig. |S| For each halo mass, M , we 
derive the probability density function of N using equation 1251 . 
Solid lines show the median occupation number as a function of 
M while dashed lines indicate the 5 and 95 percentiles of the 
distribution. 



al. 2003). It is interesting to check what are the implica- 
tions of these observational results for our halo model. We 
first note that extremely short quasar lifetimes correspond 
to very low values for Nq. From Fig.|3 we then learn that the 
condition tq > 10^ yr basically rules out all the models with 
Mq < IO'^^Mq. On the other hand, a timescale > 10^ yr is is 
fully consistent with our results for iq presented in Section 
|S] In other words, constraints to the halo model from quasar 
luminosities and from the proximity effect are consistent and 
approximately equivalent. Future determinations of quasar 
radiative histories based on the transverse proximity effect 
(e.g. Adelberger 2004 and references therein) will hopefully 
provide more stringent limits. 



6.2 Constraints from the proximity effect 

A quasar produces enhanced ionisation of H and He in its 
surroundings thus creating opacity gaps in its spectrum (or 
in the spectrum of background QSOs lying on adjacent lines 
of sight). The physical characteristics of these HII and Helll 
regions can be used to estimate the quasar lifetime (e.g. Ba- 
jtlik, Duncan & Ostriker 1988; Heap et al. 2000) . A number 
of studies have shown that, in order to explain the proxim- 
ity effect in the Ly-a forest, quasars have to maintain their 
ionizing luminosity for at least 10^ yr (e.g. Bajtlik, Duncan 
& Ostriker 1988; Schirber et al. 2004). At the same time, the 
best estimates for the quasar lifetimes indicate that tq > 10^ 
yr (Hogan et al. 1997; Anderson et al. 1999; Jakobsen et 



^^ Some authors, however, find that a disc-like component is al- 
ways needed to accurately fit the data at large radii (Percival et 
al. 2001; Hutchings et al. 2002). 



T IS THE HALO OCCUPATION NUMBER 
EVOLVING? 

The analysis presented in Section |3 is based on two ba- 
sic assumptions: the halo model and an assumed functional 
form for N{A'I), namely, equation II8II . Within these work- 
ing hypotheses, in the three-dimensional parameter space 
(q, Mo, No) we identified a one-dimensional family of models 
which accurately fits the abundance and clustering proper- 
ties of quasars in the 2QZ. Prior information on Afo, inferred 
from quasar luminosities, was used in Section 2] to remove 
the degeneracy between the model parameters. 

In this section, we want to test whether quasar clus- 
tering (without any additional constraint from quasar lumi- 
nosity) is consistent with a non-evolving model for N{AI). 
Indeed, the contours in Fig. |3 obtained for quasars at dif- 
ferent redshfits tend to lie in the same region of the a — Mo 
parameter space. This also applies to the halo occupation 
number of the entire quasar sample (0.8 < z < 2.1). We then 
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assume, as a working hypothesis, that the shape of the halo 
occupation number (parameterized by a and Mo) does not 
evolve within the time-interval spun by our quasar dataset. 
On the other hand, we let the overall normalisation A'o vary. 
In fact, because of selection effects, quasars lying at higher 
redshifts tend to be (on average) intrinsically brighter than 
their lower-redshift counterparts (cf. Table 0. Therefore, 
since we are considering objects within different luminosity 
ranges, it is reasonable to assume that they will correspond 
to different values of A'^o and, thus, to different number den- 
sities. We denote these new parameters as Nq, Nq^ and A'o^ 
respectively for the low, median and high redshift samples. 

In Fig. Q we show the confidence levels in the {a, Mo) 
plane obtained by combining the three redshift subsamples. 
The objective function (total x^) has been computed by 
adding together the x'^s of each sample. The contours shown 
in the figure are obtained by minimizing the total x^ func- 
tion over the different A^gs (we remind the reader that for 
each pair of values for (a, Mo) it is always possible to choose 
the A'oS so that to perfectly match the observed densities). 

The minimum value assumed by the total x^ func- 
tion over the parameter space is 10.98 with 10 degrees of 
freedom. Therefore, assuming Gaussian errors, our work- 
ing hypothesis that the halo occupation number of bright 
quEisars does not evolve with redshift is not rejected by 
the data at any significant confidence level. The best fit- 
ting values for the parameters are; a = 0.0^"'*, Mo = 
12.7 ±0.1, log,„iV,5' = ~1.96tl;ll logio^<f = -l-86lS:24 
and logj^Q A'^o^ — — I.TSIq'jt- ^^ ^^^ quoted intervals corre- 
spond to Ax^ ~ 1. This corresponds to 6cff ~ 2.08 ± 0.10, 
2.64 ±0.15 and 3.20 ±0.20 respectively for the low, medium 
and high redshift samples. In this case then, changes in the 
bias parameter are merely driven by the joint time evolu- 
tion of the halo population and of the mass autocorrelation 
function. 

In summary, the combined dataset is consistent with a 
model for the halo occupation number which does not evolve 
with lookback time and exhibits a very shallow dependence 
on the halo mass (a < 1, with values near zero which are 
favoured). One also finds Mo ~ 5 x IO^^Mq and A'',^ ~ 0.01- 
0.02. For all the quasar sub-samples, this corresponds to 
fQ ~ (3 - 4) X 10^ yr. 
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Figure 7. Contour levels for the x^ function obtained by com- 
bining the quasar sub-samples at different redshifts and assuming 
that the shape of the halo occupation number does not evolve 
with z. The x^ function is shown as a function of the halo model 
parameters a and Mq and it has been minimized with respect to 
A^o" , A^Q^ and A^q^ . A point indicates the best- fitting model while 
the heavy lines mark the 68.3 and 95.4 per cent confidence levels 
(respectively corresponding to Ax^ = 2.3 and 6.17). For ease of 
comparison, the contours presented in the bottom right panel of 
Fig. 121 are represented with light lines. These refer to the halo 
occupation number of our entire quasar sample in the redshift 
range 0.8 < ^ < 2.1. 



typical masses of ~ 4 x 10 Mq . The corresponding cluster- 
ing amplitude {b ~ 1.6 at z — 1.89) is too low to match our 
measures (6 = 3.9 ± 0.3 at ZcB ~ 1.89), thus suggesting that 
some revision of the model is probably required. 

On the other hand, our findings are in good agreement 
with the typical mass of haloes hosting local radio galaxies 
(Magliocchetti et al. 2004) . This further stengthens the con- 
nection between active galactic nuclei which exhibit different 
observational properties. 



8 DISCUSSION 

8.1 Comparison of results 

In Sections 14.21 13 and |7| we derived the quasar halo occu- 
pation number using a few different methods. The corre- 
sponding outcomes are fully consistent with each other. In 
all cases, we find that bright quasars are hosted by mas- 
sive haloes with M > 1O^^M0. For larger halo masses, the 
shape of the halo occupation number is not well constrained 
by the observational data and a wide range of possibilities 
is allowed. However, independently of the model details, 
we find that quasar hosts have characteristic masses of a 
few xlO^^M0. This is the key result of this analysis which 
strongly constrains quasar formation models. For instance, 
by coupling hydrodynamical simulations of galaxy forma- 
tion with simple recipes for AGN activation, Di Matteo et 
al. (2003) recently concluded that quasar hosts at z ~ 2 have 



8.2 Control of systematics 

A number of assumptions have been used in the present 
study. We discuss here how possible sources of systematic 
errors might affect our results. 

All our analysis is developed within a specified cosmo- 
logical framework based on the CDM paradigm. Modifying 
the cosmological parameters within the ranges allowed by 
recent CMB studies (e.g. Tegmark et al. 2004) induces mi- 
nor changes in our conclusions. Results similar to those pre- 
sented here are also obtained by slightly altering the power 
spectrum of density fiuctuations. For instance, neglecting 
the presence of baryons (i.e. modifying the shape parameter 
of the linear power spectrum from 0.16 to 0.21) increases 
the bias parameters of our sub-samples by ~ 7 per cent. In 
consequence, our best-fitting values for a increase by 0.2- 
0.3. At the same time, for a given a, the best-fitting values 
for logio Mo and logjQ A^o increase by 0.2-0.4. 
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The normalisation of the linear power spectrum of den- 
sity fluctuations is still very controversial: estimates of ag, 
from weak-lensing and cluster abundances range between 
0.7 and 1 (see e.g. Table 4 in Tegmark et al. 2004 for a list 
of the most recent determinations) . In Table |1] we use our 
entire sample of quasars to show how the best-fitting param- 
eters of the halo model change with irg. For simplicity, we 
just consider models with a = G and a = 1. Note that, while 
the estimated bias parameter and erg are inversely propor- 
tional, the best-fitting parameters of the halo model depend 
only slightly on the assumed value for ug (compared with 
their statistical uncertainty). 

Our analysis relies on a set of fitting functions calibrated 
against numerical simulations. These have been used to com- 
pute the mass function and bias parameter of dark matter 
haloes and the non-linear power spectrum of density fiuctu- 
ations. Considering all the uncertainties, we estimate that, 
on the scales considered here, the accuracy of the resulting 
correlation function is of the order of 10-20 per cent. This 
is still smaller than the statistical error associated with the 
observed correlation function. In consequence, we do not ex- 
pect our results to be significantly affected by this source of 
systematic errors. 

We used the most recent observational determinations 
of the Eddington ratio and of the correlation between Mbh 
and Vc to estimate the mass function of quasar host haloes. 
What is the sensitivity of our results to these assumptions? 
Assuming that all high-redshift quasars shine at the Ed- 
dington luminosity (which is a bit extreme but certainly 
plausible) would decrease our estimates for Mbh by a factor 
2-3 and the mass of the host-haloes by a factor of 3-5. The 
best-fitting solutions for N{M) would then correspond to 
values for a which are slightly larger than those presented 
in Section 2] 

A large fraction of quasar-host galaxies are morpholog- 
ically disturbed or interacting. This suggests that efficient 
black-hole fueling is triggered by galaxy encounters involv- 
ing at least one gas rich object. Based on this, Kauffmann 
& Haehnelt (2000) developed a merger-based prescription 
of AGN activation. In our analysis, we never distinguish be- 
tween merging and non-merging haloes. Can this bias our 
results? Previous studies have shown that, at z ~ 2, merg- 
ing and randomly selected haloes of the same mass have the 
same clustering properties (Kauffmann & Haehnelt 2002; 
Percival et al. 2003). This implies that our results are valid 
also in the merger-driven scenario for AGNs. However, if 
quasars are indeed found only in merging haloes, our esti- 
mates for Iq should be revised upwards by a factor of /~^r 
with /nior the fraction of merging haloes. 

8.3 Number-density evolution: implications for 
high-redshift quasars 

It is interesting to study how the number density of quasars 
with a given halo occupation number evolves. This is shown, 
for different halo models, in the left panel of Fig. |H| In all 
cases, the number density rapidly drops with redshift as 
a consequence of the hierarchical assembly of dark-matter 
haloes (see also Efstathiou & Rees 1988). On the other hand, 
by integrating the 2QZ luminosity function (Croom et al. 
2004) above a given threshold value, one finds that, be- 
tween 0.4 < z < 2.1, nqso increases with lookback time 



Table 4. Dependence of the best-fitting bias and halo-model pa- 
rameters on (Tg . The first set of data refers to models with a = 
and the second to a = 1. The entire quasar sample (0.8 < z < 2.1) 
is considered here. 
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(see Fig. IHl. This is clearly seen also in Table Q both our 
low-redshift sub-sample and our full sample roughly corre- 
spond to Mb- < —22.5 but the quasar number density at 
Zcff = 1.47 is a factor 1.3 higher than at Zeff = 1.06. We 
then conclude that, at z < 2 and for a given luminosity 
threshold, the quasar halo occupation number cannot stay 
constant with time: at least its overall normalisation, A'^o, 
(and the corresponding quasar lifetime) has to increase with 
z. This is probably due to the fast depletion of the gas avail- 
able for accretion onto supermassive black holes during the 
late stages of galaxy and group formation (e.g. Cavaliere & 
Vittorini 2000). Once again it is important to stress the dif- 
ferent nature of the halo-model parameters. Basically, while 
Mo determines which haloes are capable of hosting super- 
massive black holes, A^o and a fix the overall normalisation 
and the scaling of the halo occupation number with the halo 
mass. These two parameters are probably more infiuenced 
(with respect to Mo) by the local physics which determines 
the efficiency of gas accretion. 

It is reasonable to expect that a fixed halo occupation 
number might accurately describe the quasar density evo- 
lution at higher redshifts when gas is ubiquitously available 
within massive dark matter haloes. In fact, there is a consen- 
sus that the comoving number density of optically selected 
quasars peaks at z ~ 2 — 3 and drops rapidly at higher red- 
shifts. ^'^ This is naturally explained by CDM models where 
galaxies form at relatively late times (Efstathiou & Rees 
1988). In the right panel of Fig. |S1 we compare the density 
evolution predicted by the halo model with observational 
data from the 2QZ and the SDSS Quasar Survey (Fan et 
al. 2001). For simplicity, we set a — and we assume that 
the function N{AI) does not evolve with time. We find that 
models with Afo ^ W^^Mq are consistent with the observed 
number density of quasars with Mb < 26 in the redshift in- 
terval 2 < z < 5. On the other hand, as discussed above, a 
fixed N{M) cannot match the data a.t z < 2. Similar results 
are also obtained for a brighter sample of i-dropout objects 
detected at 2 ~ 6 (Fan et al. 2004). 

^^ A number of factors (namely, observational incompletenesses, 
uncertainties in the K-corrections and the possibility that a large 
fraction of quasars is not detectable in the optical band due to 
dust extinction) could generate a spurious drop but it is widely 
believed that at least part of the observed decrease is real (see 
e.g. the discussion in Fan et al. 2001). 
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Figure 8. Number-density evolution of optically bright quasars. Le/t panel: The solid lines are obtained using the best-fitting luminosity 
function from the 2QZ (Groom et al. 2004). From top to bottom they refer to Af^ < —22.5, —23.6, —24.1 (which correspond to the 
faintest objects in our sub-samples). The remaining lines show the evolution of nqso corresponding to a fixed halo occupation number. 
Two best-fitting models for N{M) at 2(,ft = 1-06 are represented with dashed lines: namely, the Prior B solution discussed in Section Bl 
(short-dashed) and the non-evolving model presented in Section |7| (long-dashed). The dot-dashed line shows the best-fitting (Prior B) 
model for our full sample. The dotted line is obtained by renormalising the short-dashed line so to fit the 2QZ data at 2: = 2.1 (which 
corresponds to logjg A^o = —1.65). Right panel: Datapoints with crrorbars show the high-redshift results from the SDSS Quasar Survey 
(Fan et al. 2001). The solid lines are obtained using the best-fitting luminosity function from the 2QZ (Groom et al. 2004). The dashed 
lines refer to halo occupation models of the form N{M) = A^o ■ S{M — Mq). The adopted values for A/q are indicated in the figure and 
A^O is fixed so to match the observed quasar density at 2 ^ 4. 



Knowledge of the number density evolution, however, 
does not provide enough information to determine the na- 
ture of high-redshift quasars. In fact, different halo occupa- 
tion models which are compatible with the observed evolu- 
tion of nqso correspond to wildly discrepant quasar char- 
acteristics. For instance, at z = 4, a model with a = 
and Mo = 10^^ Mq corresponds to a correlation length of 
6.7 h~^ Mpc (6 = 5.0) and to a mean host-halo mass of 
(M) = 1.8 X W^^ Mq. In this case, the observed number den- 
sity implies that iVo = 3.8x10"* andfq = 2.5x10^ yr. Thus, 
assuming that equation IIA4t still holds at z = 4, from the es- 
timated (M), one derives Mbh = 1.2 x 10* A/© which implies 
a super-Eddington accretion rate with rj = 2.8. On the other 
hand, the corresponding results for a model with a = and 
Mo = 10"Mq are: A^o = 0.18, ro = 14 h'^ Mpc {b = 9.2), 
(M) = 1.4 x W^^Mq, fq = 1.3 x 10* yr, A//bh = 2.1 x 10^ Mq 
and ri — 0.16. This clearly shows that future clustering mea- 
surements (hopefully combined with information on iq) will 
be crucial to understanding the physical properties of high- 
redshift quasars. The low abundance of optically bright ob- 
jects, however, poses enormous difficulties for this kind of 
studies. 



9 SUMMARY 

We have used a flux limited sample of ~ 14, 000 2QZ quasars 
with Mtj < —22.5 to study the quasar clustering properties 
in the redshift range 0.8 < z < 2.1. Our main results are 
summarized as follows. 



(i) For spatial separations between 1 and 20 h"^ Mpc, 
the correlation function for our whole quasar sample (cor- 
responding to an effective redshift Zcb — 1.47) is well ap- 
proximated by a power law with slope 7 — 1.53 ± 0.20 and 
comoving correlation length ro = 4.8l5'g ^~^ Mpc. 

(ii) Splitting the sample into three redshift ranges, we 
find evidence for an increase of the clustering amplitude with 
lookback time. The correlation function for quasars at 1.7 < 
z < 2.1 (zcff = 1.89) is nearly a factor of 2 higher with 
respect to the whole sample (zeff ~ 1.47). Since flux-limited 
surveys tend to select intrinsically brighter objects at higher- 
redshifts, it is not possible to tell, however, whether this 
effect is due to real evolution of the quasar population or 
to luminosity dependent clustering. We will further address 
this issue in a future paper. 

(iii) For all the sub-samples, the correlation function 
is well approximated by a power law. The best-fitting pa- 
rameters, which are strongly covariant (see Fig. 0, range 
between -2.0 < 7 < -1.5 and 4 < ro/h~^Mpc < 8 (see Ta- 
ble |5J. Within the statistical uncertainties, data in different 
redshift bins can be anyway described by the same value of 
7. Assuming that the slope of the correlation function does 
not change with redshift, evolution of the correlation length 
is detected at the 3.6 ct confidence level. 

(iv) Within the framework of concordance cosmology, 
high-redshift quasars are more biased tracers of the mass 
distribution than their low-redshift counterparts. The ob- 
served quasar-to-mass bias parameter is consistent with be- 
ing scale-independent for the separations probed by our 
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analysis. Assuming erg — 0.8, we obtain h — 2.42_g2i for 
the whole quasar sample. On the other hand, we find b = 
1.8010.24 for 0.8 < 2 < 1.3, h = 2.6215!:^ for 1.3 < z < 1.7 
and b = S.SGlg'gg for 1.7 < z < 2.1. In hierarchical mod- 
els for structure formation, the bias parameter of a popu- 
lation of tracers can be readily linked to the mass of their 
host dark-matter haloes. At a given z, values of b which 
are significantly larger than unity correspond to haloes with 
M ^ Mt{z) where Mt,{z) denotes the characteristic mass 
of haloes which are forming at that epoch out of 1 a density 
fluctuations. The bias parameters of our sub-samples then 
suggests that 2QZ quasars are hosted by rare haloes with 
M- 10^^ M0. 

(v) Using the halo model, we find that the observed 
quasar number density and clustering amplitude are con- 
sistent with a picture where: (a) quasars form in haloes 
with M > IO^'^Mq; (b) the characteristic mass of their host 
haloes is a few x 10^^ Mq. This result is independent of the 
detailed form of the halo occupation number and hence it 
can be used to constrain models of quasar formation. 

(vi) Our best-fitting models at z^h = 1.06 suggest that 
N{M) oc M°-*-°-"^ for M > W^^Mq and rapidly drops to 
zero for smaller values of M. For higher redshifts, N{M) 
tends to increase more rapidly with the halo mass. For in- 
stance, at 2cff = 1-89, N{M) oc M^'^-^ for M > lO^Mg. It 
is worth stressing, however, that the data are also consistent 
with a non-evolving functional form for N{M) where quasars 
reside in haloes more massive than 5 x 10^^ Mq and where 
the halo occupation number has a very weak dependence on 
the halo mass. 

(vii) The mean number of quasars per halo is always 
much smaller than one. Systematic searches for close pairs 
are needed to understand whether 2 active quasars can be 
hosted by a single halo. 

(viii) The observed clustering evolution is consistent 
with assuming that the locally observed correlation between 
black-hole mass and host-galaxy circular velocity (Ferrarese 
2002; Baes et al. 2003) is still vahd at 2 > 1. 

(ix) The fraction of potential host-haloes which indeed 
harbour a bright quasar increases from < 1 per cent at Zeff = 
1.06 to 5— 10 per cent at ZcB = 1.89. From this, we infer that 
the characteristic quasar lifetime iq increases with redshift 
(and/or with optical luminosity), ranging from a few xlO'^ 
yr at 2 ~ 1 to ~ 10* yr at 2 ~ 2. This is in good agreement 
with studies of the proximity effect (e.g. Jakobsen 2003 and 
references therein). 

(x) For 2 < 2, the halo occupation number of quasars 
which are above a given absolute luminosity threshold can- 
not stay constant with time. In order to match the observed 
number density evolution, at least its overall normalization 
No (and thus the corresponding iq) has to increase with 2. 
This probably refiects the fast depletion of the gas avail- 
able for accretion onto supermassive black holes during the 
process of galaxy formation. 

In brief, this paper presents state-of-the-art measure- 
ments of quasar clustering and establishes an accurate 
benchmark for quasar-formation models. Future results from 
the SDSS quasar survey will provide an independent veri- 
fication of our results and, thanks to the different quasar- 
selection criteria, will extend them to even higher redshifts. 
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APPENDIX A: THE MASS DISTRIBUTION OF 
QUASAR HOST HALOES 

In this section, we use a few observational results to derive 
the conditional probabiliy distribution of the host-halo mass, 
M, for a quasar with given absolute magnitude M^j . 



Al From photographic to Johnson B magnitudes 

We start from converting 6j fluxes into standard B magni- 
tudes. In general, B ~ 6j + 0.3 (B — V^) (Blair & Gilmore 
1982; CoUess et al. 2001), and the rest-frame color index for 
quasars is B - V ~ 0.22 (Cristiani & Vio 1990). In what 
follows, we then assume 



Mb 



Aft, 



0.07 



(Al) 



which is in good agreement with Brotherton et al. (2001). 
Note that the amplitude of the correction is comparable with 
the statistical error which affects the magnitude determi- 
nation in the 2dFGRS (CoUess et al. 2001; Norberg et al. 
2002b) which uses the same UKST photographic plates as 
the 2QZ. It is then reasonable to expect that also quasar 
photometry in the 2QZ is affected by typical errors of ~ 0.1 
magnitudes (e.g. Corbett et al. 2003). 



A2 Bolometric corrections 

In order to use equation 1201 1 to infer the mass of the black 
holes which power 2QZ quasars, we need to convert their 
absolute B magnitudes into bolometric luminosities. Bolo- 
metric corrections for a sample of X-ray selected quasars 
lying at 2: < 1 have been derived in a seminal paper by Elvis 
et al. (1994). Since observations show that quasar spectra 
do not evolve with z (e.g. Bechtold et al. 1994), it is com- 
mon practice to apply these corrections also to high-redshift 
quasars. It has been recently pointed out, however, that the 
bolometric corrections by Elvis et al. (1994) are seriously af- 
fected by systematics and should be revised downwards (e.g. 
Fabian & Iwasawa 1999; Elvis, Risaliti & Zamorani 2002). 
For this reason we use here the results by McLure & Dun- 
lop (2004) who, adopting the revised template spectrum by 
Elvis et al. (2002), estimated the bolometric corrections for 
1136 quasars at 0.5 < z < 0.8 extracted from the SDSS. 
Corrections from the B band have then been computed us- 
ing a subsample of 372 objects common to the the 2dF and 
SDSS surveys. When combined with equation lAlL their 
best fltting relation gives 



( lo4rgs-0 ^ °-'' " °-'' ^^'''" + ''^ ^"^'^ 



and the corresponding rms variation at fixed Mb^ is 0.14. 
Within the quoted uncertainties, this is perfectly consistent 
with the recent results by Marconi et al. (2004). ^* 



A3 The Eddington ratio and the distribution of 
black-hole masses 

Observational estimates of the Eddington ratio, rj, require: 
(i) using some dynamical tracer to determine the black-hole 
mass (and thus the Eddington luminosity); (ii) measuring 
the quasar luminosity in a given band, Li; (iii) applying the 
corresponding bolometric correction, /3f° ; (iv) calculating 
77 (X /jj'"' Li/Mbh. Given this complexity, measurements of 
rj are rather uncertain and sensitive to a number of sources 
of systematic errors. Recent determinations, however, tend 
to lie in the same ballpark and suggest that r] mildly in- 
creases with z (Dunlop et al. 2003; McLure & Dunlop 2004). 
For consistency with Section rA2l we use here the results by 
McLure & Dunlop (2004) who combined virial estimates of 
black hole masses with new bolometric corrections to infer 
the Eddington ratio for a large sample of quasars from the 
SDSS. From their results we infer that, for 0.8 < z < 2.1, the 
probability density function for logj^Q ry is well approximated 
by a Gaussian distribution with mean 0.21^ — 0.80 and vari- 
ance ~ 0.3. This corresponds to {7?) = iq0.21z~q.65 ^ These 
results are also supported by other indirect determinations 
of rj. By requiring the mass function of relic black holes (as 
inferred from the X-ray background) to match its local coun- 
terpart, Marconi et al. (2004) found that 0.1 < r? < 1.7 (with 
a preferred value of 7; ~ 0.5). Once accounted for the dif- 
ferent bolometric corrections, these values are in extremely 
good agreement with the results from McLure & Dunlop 
(2004). Similarly, Yu & Tremaine (2002) showed that the 
local mass density in black holes is consistent with the inte- 
grated luminosity density of quasars if they accreted mass 
nearly at the Eddington rate at redshifts 2 > 2. 



A4 The distribution of black-hole masses 

The conditional probability distribution of logj^Q Afbh/M© 
for a given Lboi is thus obtained by combining equation 1201 
with the observationally determined distribution of 77: 



^(logi 



Mb, 

M,T 



\Lhoi) = 0.73 exp 



(logio^-/) 



0.6 



(A3) 



with / = 8.70 - 0.212 -h logio(iboi/10*'^ergs"^). For con- 
sistency, in order to estimate the black-hole mass associ- 
ated with a quasar of given absolute magnitude Mf,j, we 
combine equations l|A2fl (including its associated scatter) 
and llA3l l which have been derived from the same dataset. 
This implies that a quasar with A/tj = —25 corresponds 
to a mean black-hole mass of 5.13 x 10* M0 at 2 = 1 
and of 3.16 x 10* M© at 2 = 2. Note that for the typ- 
ical redshift and magnitude ranges spun by our sample, 
6 X 10^ Mq < (MbhlMtj) < 3 X 10^ Mq. This interval is con- 
sistent with the masses inferred from dynamical measures in 
the local Universe (e.g. Tremaine et al. 2002 and references 
therein) and from the analysis of emission linewidths in the 
2QZ (Corbett et al. 2003). 



^^ The bolometric corrections by Marconi et al. (2004) are 
roughly 2/3 of those by Elvis et al. (1994) and correspond to 
logio ( io-ibV°'s-0 = ^■^'^ ~ 0.40 (Mtj -I- 25) with a scatter at 
fixed Mb, of 0.3. 



A5 The mass of host haloes 

Taking a step further, we can estimate the probability dis- 
tribution that a quasar of a given luminosity is hosted by a 
dark-matter halo of mass M. Following Ferrarese (2002; see 
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also Baes et al. 2003), this is obtained by assuming that a 
statistically significant correlation links Mbh and M. 

Black-hole masses are found to be tightly correlated 
with the velocity dispersion of their host spheroid, (Tsph (Fer- 
rarese & Merritt 2000; Gebhardt et al. 2000). The most re- 
cent determination considers ~ 30 galaxies with secure de- 
tections of supermassive black holes (Tremaine et al. 2002) . 
Observations also provide evidence for a correlation between 
(Tsph and the circular velocity in the flat part of the rotation 
curve of the host galaxy (Ferrarese 2002; Baes et al. 2003). 
By combining the Afbh — o"sph and the asph — Vc relations, 
Baes et al. (2003) find 






10 



7.24±0.17 



200 km s- 



(A4) 



This purely observational relation can be used to link Mbh 
with the mass of the host halo. In order to do this, however, 
one needs to express Vc in terms of M which is a formidable 
task. As a first order approximation one can assume an equi- 
librium configuration for the dark-matter density profile in 
haloes. Both the singular isothermal sphere and models de- 
rived from numerical simulations (e.g. Navarro et al. 1997) 
provide good starting points. However, detailed modelling 
of the rotation curve requires accounting for the distribu- 
tion and physics of baryons (e.g. Mo, Mao & White 1998). 
In fact, the gas contribution can be dominant in the in- 
nermost regions of galaxies. Moreover, the condensation to- 
wards the centre of the dissipative material redistributes, 
through gravity, the coUisionless matter. 

For simplicity, we consider equilibrium profiles which 
only contain dark matter and we account for the presence 
of baryons in an approximate way. The circular velocity at 
the virial radius of each halo is 



159.4 km s-i 



M 



1012 /i-i Me 



1/3 



18 7r2 



1/6 



(A5) 



where E'^ = Qq (^ + z)'^ + f^A, A^ is the ratio between the 
mean density of the halo and the critical density of the Uni- 
verse (both evaluated at redshift z). For a spherical collapse, 
this function can be approximated as A ~ 18it^ +82x ~S9x'^ 
with X = f7m(z) - 1 and ^^(z) = ^o (1 + zf jEl (Bryan & 
Norman 1998). 

A truncated singular isothermal sphere has a constant 
circular velocity profile v^ = t^vir, while for an NFW density 
profile 



ATI) 



1 F{C ■ 7^) 



7^ F{C) 



1/2 



(A6) 



where Ti. — r/rvh, F{x) = ln(l -\- x) — x/{l + x) and C is 
the concentration parameter of the halo. In this case, the 
circular velocity vanishes when 7^ ^ 0, reaches a maximum 
at 7?. ~ 2.16/C and matches the virial velocity &tTZ = 1. Two 
questions naturally arise: t) What is the value of IZ which 
corresponds to the observed circular velocities? li) What is 
the contribution of the baryons at this radius? These are the 
main uncertainties of our analysis. 

For galaxies with HI rotation curves, Vc is typically mea- 
sured at a few tens of kpc from the centre, well beyond the 
optical radius (a few kpc). On the other hand, the present- 
day virial radius of a halo with M = lO^^M© is rvir = 0.56 
Mpc. In other words, the largest scales sampled by rotation- 
curve measurements are nearly a factor of 10 smaller than 



the virial radius. Using galaxy-galaxy lensing data from the 
SDSS, Seljak (2002) has shown that, for galaxies above L,, 
Vc decreases significantly from the optical radius of a galaxy 
to the virial radius of its host halo. This result is indepen- 
dent of the morphological type and is probably suggesting 
that baryons contribute significantly to the circular velocity 
at the optical radius and that density profiles for the dark 
matter are highly concentrated (as expected in CDM models 
a,t z = Q). Seljak (2002) also found a clear trend for the ra- 
tio i/j = Vc/v^ir with halo mass. Typical values are: tp ~ 1.8 
for A/ ~ 3 X 10" Mq, ^ ~1A± 0.2 for M ~ IO^Mq and 
-i/j < 1 for cluster masses. This is in good agreement with the 
predictions of CDM models, since the dark matter concen- 
tration is expected to decrease with the halo mass and the 
baryonic contribution is expected to become less and less 
important. 

Assuming that the observed Vc corresponds to the max- 
imum value of the rotational velocity profile in an NFW 
halo tends to underestimate Seljak's results. Using equa- 
tions (9) and (13) in Bullock et al. (2001), we find that, at 
2 = 0, this assumption corresponds to i/; = 1.3 (C ~ 14) for 
A/ ~ 3 X 10" A/0 and ip = 1.2 (C ~ 9) for M ~ IO^^A/q. 
Anyway, these results show the correct trend with the halo 
mass: smaller, more concentrated haloes are associated with 
larger values for tp. It is worth noticing, however, that the 
candidate hosts of our quasars (haloes with M ~ W^^Mq 
at 0.8 < < 2.1) are expected to be much less concentrated 
(C ~ 3 — 5) than their present-day counterparts. In this case, 
the maximum value of the rotational velocity is only 2-15 
per cent higher than iivir- This motivates the choice -(/' ~ 1 
as a viable alternative to the low-redshift results by Seljak 
(2002). 

We have now collected all the elements necessary to es- 
timate the conditional probability distribution of the host- 
halo mass, M, for a quasar with given absolute magnitude 
Mb , : P (M I Mfc , ) . In brief: i) We assume that equations 1A2II , 
<A3l l and llA4l l. which have been determined at lower red- 
shifts, are still valid for the host galaxies of our 2QZ quasars 
at 0.8 < 2 < 2.1; Their combination (including the scatter 
in each of them) is used to determine the probability distri- 
bution oivc; a) We then convert circular velocities into halo 
masses by selecting a value of ip and using equation <A5t . To 
match the results by Seljak (2002), we assume a Gaussian 
distribution for tp with mean value of 1.4 and scatter of 0.2 
(case A); alternatively, based on the estimated low concen- 
tration of high-redshift haloes, we assume that tp — 1 (case 
B). 



APPENDIX B: THE HALO OCCUPATION 
DISTRIBUTION FROM SEMI-ANALYTIC 
MODELS OF GALAXY FORMATION 

We use here semi-analytic models of galaxy formation to get 
an insight into the problem of choosing a functional form 
for the first two moments of the quasar halo occupation 
distribution. 



Bl The halo occupation number 

There is evidence that, at high redshift, quasars are asso- 
ciated with star-forming galaxies (e.g. Omont et al. 2001; 
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Hutchings et al. 2002). It is then plausible to expect that 
the halo-occupation properties of quasars might share some 
similarities with those of galaxies which show active star 
formation in their nuclear region. As an example, we de- 
rive here the function N{M) from the semi-analytic mod- 
els of the GallCS I collaboration (Hatton et al. 2003) at 
z = 1.08 (roughly corresponding to the median value for our 
low-redshift sample). Results for galaxies with a bulge star 
formation rate, i/Jbuige, which is larger than lOM0yr~^ are 
shown in Fig. IBll The choice of such a threshold for V'buigc 
is motivated by fact that the mean density of these objects 
(~ 12 X 10~® h^ Mpc~^) is comparable with the mean den- 
sity of our low-redshift quasar sample. In order to improve 
the statistics, in Fig. lBll we also show the function N{M) for 
galaxies with i/Jbuigc > 2 Mq yr"'^. In both cases, the halo oc- 
cupation number is well approximated by a power law with 
a cutoff at small virial masses. For instance, the function 



N{M) = No X 




if M > Mo 



(Bl) 



very closely matches the results of the semi-analytical mod- 
els in Fig. IBll This is in good agreement with equation 
11811 where a sharp cutoff replaces the exponential decline 
at small masses. 



B2 The scatter of Pn{M) 

In this section, we use the previously introduced samples 
of star-forming galaxies to study the second moment of the 
halo occupation distribution. We first note that there is not 
a single halo in the GallCS I sample which contains more 
than 1 galaxy with iphuigc > lOM0yr~^ at z = 1.08. In 
other words, the data (within extremely large errorbars) 
are consistent with E^ = 0. This is clearly an effect of the 
small number statistics. On the other hand, the results for 
V'buigG > 2M0yr~^ are much more significant. In this case, 
the numerical results are well approximated by the function 



r(M) 



ii 
Ms 



1 + 



Me 



(B2) 



which scales as a power law for M <C M^ and approaches 
1 for M ^ Mb (see Fig. IB2II . In agreement with studies of 
low-redshift galaxies (e.g. Sheth & Diaferio 2001; Berlind & 
Weinberg 2002; Berlind et al. 2003), the scatter of Pn{M) 
is then strongly sub-Poissonian for haloes which, on aver- 
age, contain less than 1 object and nearly Poissonian for 
larger haloes. We find that, adopting N{Ms) — 0.75 as an 
operative definition for Ms, equation <B2ll with 7s = 2 accu- 
rately describes the second moment of the halo occupation 
distribution for rare galaxies at z ~ 1. This result does not 
depend on the details of the galaxy population considered 
(star formation rate, colour, etc.). 



-2 - 




13 

log,„(M/Mj 

Figure Bl. The halo occupation number of galaxies which at 
z = 1.08 are actively forming stars in the bulge as obtained from 
semianalytic models of the GallCS collaboration. Datapoints cor- 
respond to the estimated N{M), vertical errorbars mark the as- 
sociated 1 a uncertainties (assuming Poisson statistics for both 
the number of galaxies and the number of haloes in a bin) , while 
the horizontal errorbars denote the size of the mass bins. Arrows 
mark the upper limit for N(M) in the bins where we measure 
N = 0. The continuous lines show a fit to the data obtained by 
using the function in equation IBll . For Vbulgo > lOAf0yr~^, 
the best-fitting parameters are a = 0.85, Mq = 10^^ M0 and 



A^o^ 



X 10" 



On the other hand, for Vbulge > 2 Mq yr , one 



gets a = 0.5, A/q = W^^'^^Mq and No = 0.42. 




Z=1.08 

V'buige>2 Mg/yr 
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Figure B2. As in Fig.|^but for the function r(Af). The fitting 
function represented with a continuous line is given in equation 
IB2i and corresponds to the best-fitting parameters 7s = 2 and 
A/s = 10^^'^Mq. The shaded region indicates the mass range 
where the function F is totally undetermined. 



